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Abstract. We introduce a weak formulation for a system of electrostatic and hydrodynamic 
equations modelling the electrophoretic motion of charged particles in ionized fluids. We obtain 
a local in time existence theorem, using the results established in [11] and properties of the 
solutions of the Poisson-Boltzmann equation. These properties follows from singular integral 
operators techniques. 

Mathematical subject classification: 35J60, 92C05, 35J50, 42B20. 

Key words: Poisson-Boltzmann equation, singular integral operators, Stokes problem, elec- 
trophoresis. 



1 Introduction 

In this paper we establish results on existence and uniqueness of the weak solu- 
tions for a system modelling the motion of a charged particle driven by the action 
of an external electrical field. This phenomena is known as electrophoresis and 
is important in many technical applications (a vast literature is available: see for 
example [29], [24], [27], [2], [1], [3], [13], [30], [14]). 

We are considering a particle (a charged polymer, for example) immersed in 
an ionized solution (a viscous incompressible fluid). On the boundary of the 
enclosure an external electric field induces the electrical potential inside the 
enclosure which is determined by the Poisson-Boltzmann equation (see [19], 
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2 WEAK SOLUTIONS FOR THE ELECTROPHORETIC MOTION 

[10]). The hydrodynamic behavior of the system is governed by the Navier- 
Stokes equation. 

We are not considering the usual approximation for the effect of the electrical 
field on the particle based on Prandtl boundary layer theory well known in colloid 
science, the so called slip-velocity condition (for details see [4], [22], [27]); their 
derivation requires better regularity properties of the boundary of the particle 
(see [28] and the discussion in the introduction of [27]). 

As remarked in [24] the theoretical analysis of the electrophoretic motions is 
quite difficult, as it combines specific features of polymer physics with the intrin- 
sic complexity of electrokinetic phenomena. From a mathematical standpoint 
the difficulties reside in the treatment of the electrical-hydrodynamic couple and 
on the low regularity of the data. We have established elsewhere [6] the existence 
of a //'-variational solution for the electrostatic potential, for a general class of 
domains. In the case of Lipschitz regions we have established a H 3/2 -regularity 
result by means of the singular integral operators theory; this regularity is opti- 
mal, even in C 1 -domains (see the comments and negative results in [20]). For 
C 1 "-domains, 0 < a < 1, and suitable charge distribution this theory can be 
applied in the classical sense ([25], [18], [7]) in order to obtain more regularity 
for the potential [5]. 

Recently, the motion of rigid bodies in a bounded domain filled with a viscous 
flow has been treated rigorously ([11], [17]). Special techniques (from the trans- 
port theory [23]) have been used in order to obtain existence and properties of 
the suitable weak solutions for these systems. In particular in [1 1] a global weak 
formulation is introduced and existence of solutions local in time is established 
when a L 2 body force and C 1,1 -domains are considered. Evidently, in the study 
of the electrophoretic motion we can not use directly these results because we 
have an external electrical field interacting with the ionized solution. However, 
we obtain a similar result of local existence (see Theorem 3.1); this is obtained 
as a consequence of uniform bounds and convergence properties involving the 
electrical force term F (Corollary 4.1, Theorem 4.2). Following the discussion 
in [6] we choose to prove these properties restricted to the case which the surface 
charge distribution of the particle and the fixed charge distribution (in the parti- 
cle and in the solution) are L 2 and L°° functions, respectively; in this case we 
need only consider the Lipschitz regularity on the boundary of the particle. The 
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existence of weak solutions obtained in Theorem 3. 1 follows from the additional 
hypothesis that the particle domain is of class C 11 . More properties on F can be 
obtained in the C l of context if we assume stronger hypotheses [5] on the charge 
distributions but we do not consider this situation in this paper. 

2 The governing equations 

Consider a rigid, charged particle immersed in a electrolyte solution under the 
action of an external electrical field. We suppose that the solution (a viscous 
fluid) occupies a region DcK 3 and at the initial moment of time, the particle 
(a rigid body) occupies a compact region Kq c D such that its center of mass is 
located at the origin y = 0 of a Cartesian coordinates y. 

Let us define K(t) as the domain occupied by the particle in the time t and 
0(x, t) as a real valued function which represents the electrical potential in 
(x, /) 6 D x [0, 7], where V/ e [0, T], distOD; dK(t)) > d, where d is a 
fixed constant. We set \J/(x, t) = ( * >( ^' )< ' , where A is the temperature (constant) of 
the system and e is the electron charge. The governing equations and boundary 
conditions to \f/ are (see [6], [19]) 

V • (k(x, t)VxJ/(x, /)) - b(x, xj/(x 7 /)) = p(x, /), x e D, 
fc(x, ')-*<*.'). x€6tf(f), 

to(x,f) = ¥(x), xedD, (21) 

d\J/\ d\I/2 4;r e 

*i-^(x, t) - fe-fV t) = — <r(x, 0, x e dK{t). 
3/1 dn A 

Here 

• k:Dx[0J]-> L(R 3 , R 3 ) is defined as k tj {x, t) = <$„*, if x € K(t), 
k u (x, t) = 8ijk 2 if x € D\K(t), where k u k 2 are the dielectric constants 
in K (t) and D\K(t) respectively. 

• b : D x R -+ R, b(x, \J/(x, t)) = k 2 r~ 2 sinh if/(x, t) if x € D\K(t), 
b(x, V(x, f)) = 0 if x € K(t), r^ 2 is the Debye radius [10]. 

• *(x) = Vn(x, 0 = i^(x, 0|jr ( », and ^ 2 (x, /) = V(x, t)\ DxW - y 
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• a is a superficial charge distribution and p(\, t) = (p\(x, t), p 2 (x, /)), 



are the charge distribution in K(t) and D\K(t) respectively. 

The action of the electrical field on the particle produce its motion. The motion 
of fluid is described by the velocity field v / (x, /) (velocity of the fluid material 
point which has Cartesian coordinates x at time t) and satisfies the Navier-Stokes 
equation 

v f (d t \ f + div(v' ® v / )) - //Av / + Vp = yjrF, in V(Q T ) 3 
div v' = 0, in Qt 

(2.2) 

\ J = 0, in 3D 

v / |, =0 = v 0 / in D\K 0 

for all / 6 (0, T). Here t] > 0 is the viscosity of the fluid, v/ > 0 is the 
homogeneous fluid density (of the mass) and Q T — {(', x)/' € (0, T), x € 
D\K(t)}\ denoting p' ons as the ion density of the solution, we have F = -(p% + 
p ions )V<f> 2 as the electrical force on the fluid domain (see [30]). Then F = 
4^ (P2 + r £> 2 ^2 sinh (\J/ 2 )) (V^2)/d\?(0' us * n 8 tne Boltzmann distribution for 
p ions (see [16]). 

Let us set x c (/) as the center of mass of the particle; w(f ) the rotation vector; 
R(0 the translational velocity; A the symmetric inertial matrix; \ p the velocity 
of the particle. We observe that if v p > 0 is the density (of the mass) of the 
particle, 



for ally € R 3 . We have also \ p (\, t) = R(f) + w(f) x (x-x c (/)) forx € K{t). 
It is important to observe the implicit dependence of F on v p . 

From the Newtonian mechanics for rigid bodies and the stress tensor in fluid 
dynamics, if Af is the mass of the particle, the evolution law for the motion is 
given by 



where 



P,(X, /) = - — P?(X), 02(X, 0 = -^-/°2°( X > 0. 

p°\ = Pi = p°\ D \m 




M 



clR(t) 
dt 



JdK(t) 



f or*(x, 0 • n(x, t)ds + 



BKU) 



a fc (x,r) n(T, t)ds 
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and 

A dMo = f (x _ xAt)) x (flr « (Xt r) . n(x f))</J + w(0 x {A , w)(0 

+ f (x-x c (0) x (a £ (x,r)-n(x,0)^. 

If we set D(v^) = i(Vv / + (Vv^) 7 ') thena^x, f) = 2r?D(v^(x, r)) - p(x, r)I 
is the stress tensor of the fluid; 

is the electrostatic tensor (see [30]). 
We assume the following hypotheses on the data 

(i) tf 0 is a Lipschitz domain. 

(U) DisaC 2 -domain,vI/ g //'(aD)nCOD),cr(. , /) € £, 2 (dJT(*)), p(. , *) € 
L°°(D),V/ e [0, 7]. 

Remark 2.1. As remarked in the introduction of this paper, under the above 
hypotheses we have established elsewhere [6] that 

tfr(. , f) = (xf; u xfr 2 )(. , f) € H\D) fl (H 3/2 (K(t)), H 3/2 (D\KiF))), 

for all / e [0, 7]. This regularity result and Trudinger's inequality (see discussion 
in [6]) give us that 

|| sinh (V 2 )(. , t)\\o. P .D\m < +°°< v 1 < P < +°° • 

Recalling the Sobolev embedding H ]/2 (D) C L 3 (D), the Holder's inequality 
gives us 

/ sinh 2 (V 2 (x, 0)1 Vir 2 (x.t)\ 2 dx 

JD\K{t) 

< || sinh (* 2 (. , OM^^I^C , 0llJ t3>oxTO < +oo. 

Then F(. , t) e L 2 {D)\ V/ € [0, T] (evidently this implies F € L 2 ((0, T) x 
D) 3 ). AsimilarcalculationshowsusthataJ(. ,r) € L 3/2 (D\K(t)),Vt 6 [0, 7]. 

As can be seen in the following section an existence result local in time of 
the suitable weak solutions for (2.2) coupled with (2.1) is available if we as- 
sume that K 0 is a C 11 -domain. 
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3 The notion of weak solution 

Following [11], let us define the Eulerian densities v p (x,t) = v p I K{t) (x), 

v f (x,t) = ^// D \^(7)(x) and the global density v = v p + v f . We also define 
the global velocity in D as 



u(x, 0 = 



?/(x,0 if xeD\K(t), 
v'(x,/) if xeK(t). 



In view of the conservation of mass, v satisfies the linear transport equation 
inD 

B,v + div(vu) = 0. 

We require that W ■ n = yf • n and o H • n = T in dK(t), where -T is the 
force applied by the particle on the fluid. We can write T = £ n, where £ is 
the Cauchy stress tensor in the body. 

On the walls, we enforce homogeneous Dirichlet boundary conditions u\ap = 
0. Moreover, the incompressibility of the fluid, the rigidity of the structure and 
yP ■ n = yf • n imply that div u = 0. 

The evolution laws of the momentum for the fluid and for the particle are 
given by 

d,(v f u) + div(v / u<g>u) = ^-div(y/(2^D(u) - pi)) + • Vv- + v/F 

1 1 U 1 p 

d,(v p u) + di\(v p u (g» u) = — div(v p i;) - —o H ■ Vv p - —o E • Vv p , 

v p V p Vp 

respectively. Here D(u) = \ (Vu + (Vu) r ) is the global rate-of-deformation 
tensor. 

Introducing the global stress tensor 



// 



J — + 



VpY, 



v f v p 

we obtain the global system in D'((0, T) x D) 3 , 

8,(yu) + div(vu <8> u) = div T + vF, 

(3.1) 

div u = 0, d, v + div(vu) = 0, v p D(u) = 0, 
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a(x,0)H V ° /(X) ' XeD ^° (3.2) 
v*(x, 0) = w(0) x (x - x c (0)) + R(0), x € K 0 , 

v p (x, 0) = Vp/ir 0 (x), v,(x, 0) = V// D ^(x) 

v(x, 0) = !> 0 (x) = v p (x, 0) + v f (x, 0). 

We reproduce here, following the paper [1 1], the notion of the weak solution 
of the above system 



Definition 3.1. (v, u) is a weak solution of (3.1H3.3) in (0, T) if it satisfies 
the a priori energy bounds 

v e L°°((0, T) x D), u € L°°(0, T; L 2 (D)) 3 n L 2 (0, 7; H*(D))\ 

and if for all <f> e V and for almost every r g (0, !T), 

(vu • 9,<p + vu <g> u : D(^) - r)D(u) : D(<p) + vF • (p) dxdr 

+ / v 0 Uo • <p(0)dx = I / vu (pdx J (/), 

• /d VJd / (34) 

9,v + div(vu) = 0, div u = 0, 

v p D(u) = 0, u| aD = 0, in T>'((0, T) x D) 3 , 

v 0 e L°°(D), u(.,0) € L 2 (D)\ 

where V is defined by 

y={(pe H\(0, T) x D) 3 /cp(t) e V(t), Vf € (0, T)} , 

and 

V(0 = {cpe // 0 '(D) 3 /div <p = 0, v p D(cp) = 0} . 
The following existence theorem for the above weak solutions is available [11] 



/7 

JO JD 



Theorem 3.1. Under the hypothesis (ii)-(v) (see Section 4) and the additional 
assumptions that Kq is a C 1,1 -domain, Uo e Hq(D) 3 , div uo = 0, v p D(uo) = 0 
and 8(0) > d, there exist T* € (0, +oo] and a solution (v, u) of (3.4) such that 

(i) 0(v) € C([0, r]; n L°°((0, oo) x D) for all T < T*, p < oo 

and jSeCHK; R). 
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(ii) u € L°°(0, T\ H*(D)) 3 and d,u e L 2 ((0, T) x D) 3 for all T < T*. 

In [11] the hypothesis F e L 2 ((0, T) x D) 3 (a body force) is assumed (see 
Remark 2. 1) and the proof of the existence theorem as Theorem 3. 1 is based on 
the solution of an approximated system (obtained by regularization techniques). 
The existence of the approximated solutions is obtained by the Schauder fixed- 
point theorem (see [5]), via a solution of an appropriate inhomogeneous (linear) 
Stokes equation. Using the C 11 -regularity of the domains and the smoothness of 
the coefficients, this linear problem has a solution with the necessary regularity. 
The solution (u, v) is built as a limit of these approximations; the existence of this 
limit is derived from the compactness properties of the linear transport equation 
[12]. This is made possible if we can obtain elliptic estimates and a priori 
bounds for u as well as energy bounds for v (see Section 4 in [11]). However 
as F depends on u we need to take some care in this regard. More precisely 
let us define for each m e N, (u {m \ v {m) ,F (m) ) such that v (m) is bounded in 
L°°((0, r)xD)uniformlyinm,u (m) is bounded inL 2 (0, T; H^(D)DW lA (D)) 3 , 
d t u (m) is bounded in L 2 ((0, T) x D) 3 uniformly in m, v { 0 m) converges to v 0 
in L 2 (D), u ( 0 m) converges to u 0 in L 2 (D) 3 and (3.4) is valid, for all <p (m) € 
y (m) . We set K (m) (t) = M im) (t)K 0 , where M (m) {t) is an invertible affine 
transformation, and we suppose that 8 = inf{6 (m) (r), t e [0, T], m > 0} > d< 
where 8 {m) (t) = distOD, dK {m) (t)); F (m) is defined by the calculation of x// (m \ 
the solution of (2.1) considering (K (m) (t), p (m \ o (m) ). We have to show that 
f 0 ||F (m) (. , t)||q 2 D dx < C, where C does not depends on m. 

Admitting this uniform bound on F (m) , the stability results [12] for linear 
transport equation can be used as in [11]: there exist (v, u) such that up to the 
extraction of a subsequence, p(v {m) ) converges to fi(v) weak* in L°°((0, T) x D) 
andinCQO, T]\ L P (D)) for all p < +ooandall0 e C 1 ^), and u (m) converges 
to u in C([0, T]; //J(D)) 3 for all s < 1. However, if F is related with \fr, where 
\fr is the solution of (2.1) considering K(t), p, cr, it is not obvious that (3.4) holds 
for (u, v, F) for all given <p € V. This is established by the special argument in 
Section 4 of [1 1] if we can show that 

f f ¥ m • (pdxdx -> / f F • (pdxdz, m -+ oo, V<p e V, Vr € [0, T] . 
Jo Jd Jo Jd 

As we shall see in the next section, in order to establish these results for F (WI) , F, 

we need to study the properties of the solution of (2. 1). Additional hypotheses 
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on p and o are also necessary and as remarked in the Section 2 we shall prove 
the results in a more general framework: we assume hypothesis (i), i.e., Kq is a 
Lipschitz domain. 

4 Bounds and convergence for the potentials 

As remarked in [6], the determination of the charge densities in bio-molecular 
systems is a non-trivial question which is treated in some computational studies 
using Hartree-Fock approximation techniques (see [9], [26]). Here we assume 
explicitly the following hypothesis 

(iii) Let p (m) as in the Section 3. Then \\p {m) \\v*mT)*D) < C, where C does 
not depend on m \ 

(iv) For all t e [0, T], o (m) {M (m) {t)x, t) = a(M (m) (t)x, t) = a(M(t)x, t) = 
<t(x, 0),VxedK 0 . 

(v) ||p (m) - p||ico ((0t DxD) 0, m -► +oo. 

Recalling that V e H l (dD) and following [6] we consider the problem 

V • (k(x, /)V$(x, /)) = 0, x e D 
$(x, t) = * (x), x € 3D, 

which has a solution $ g H 3 ' 2 (D) such that %-,^e L 2 (dK(t)), V/ e [0, T] 
and _ 

^(•,/) < C||$(.,Oll3/2.2.D, (4.1) 

II dn II 0.2.cl/C(r) 

where C depends only on the Lipschitz nature of dKo (see the papers [8], [31] 
for details). We observe also that $(. , /) g L°°(D) (see [15] or [21]). 
Introducing V = *A + we see that (2.1) may be reformulated as 

V • (k(x. f)V?(X, /)) - Mx, ?(X, f) + $(x)) = p(x, /). X G D, 

fc(x,0 = ^i(x.r). xedK(t), 

fe(x, /) = 0, x e 9D, (4 2) 

d\j/\ d\jr~> 4ne ( 3vpi 3$-> \ 

k l ffix,t)-k 2 -jj2(x,t) = _-a(x,0-^,^-(x,0-*2^(x,0j 

= a(x, /). on aAf(0- 
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Let us consider t e [0, T] fixed and recall the weak formulation for (4.2). 
Find ueV = H* (D) such that b(x, u + $) € L 2 (D) and satisfying 

a(u, v) + (#(«), u) = L(w), Vti 6 , (4.3) 

where 

fl(n,i;) = f k(x,t)VuVvdx,(N(u),v) 

JD 

= / b(x, u + $)vdx, L(v) 

JD 

= / ffftwfo - / /ow/x, Ko 

•WO JD 

is the usual trace operator. Here k(x, t ) = k\ if x € K"(0> £(*» 0 = ^2 if 
x € D\KX0. 

As can be seen in [6] (see also [19]), this problem is equivalent to find the 
minimum in Hq(D) of the functional 

F(u) = \f 0\Vu\ 2 dx - L(u) + J(u) (4.4) 

where /(.) is defined as 

J{u) = k 2 rp 2 I {cosh (k + $) - cosh ($)}dx, 

JD\K(t) 

if / | cosh (m + $) -cosh($)| 2 dx < oo, 

JD\KU) 

J(u) = + oo if the square integral is + oo. 
Below we establish the first bound on \f/ derived from (4.3) 

Lemma 4.1. Let us assume the hypothesis (i)— (iv). Then the solution \j/ e 
H\D) 0/(2.1) belongs to L°°(0, T; H { (D)) and UWl^oj-hHd)) < M, 
where M depends only on ||/o||l~((OT)xD), , f)Ho,2,dr<f)» k h k 2 , r~ 2 , D, 
dK 0 andW. 
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Proof. By \f/ = \J/ + 4> we only need to prove the lemma for Using (4.3) 
we have, for all t € [0, 7*], 



-*] f |V#(x, t)\ 2 dx -k 2 [ _ |V^(x, 0| 2 ^x + f (o>o£)(x, /)<fc 
./ku) Jn\Ku) Ji)K(t) 

= f (/oi?)(x,/Mx + ^ 2 r- 2 f sinh(£(x,/) + $(x,0W(x,/)dx 

or 

min(*i,*2) /" |V^(x,/)| 2 </x 

< / (ayoh(xj)ds - f (p$)(x,t)dx+ 
JnK{i) 3d 

- k 2 r- D 2 f sinh(£(x,r) + *(x. 0)(£(x. 0 + $(x./))</x 

+ *2^ 2 / _«nh(^(x./) + 4mx.OWx,/)Jx 

< IIP(.,0llg, 2 , D «ill?(..0llg, 2 , f) ll^-^)ll( 2 ), 2 .^(n € 2\\(Yoh(.,t)\\l 2 , i)KU) 
+ k 2 r- 2 \\$(..t)\\ 00 . D [ |sinh(£(x,/) + Y(x./))|</x. 

where we have used inequalities of Schwarz and Young. Now, recalling that \j/ is 
the minimum of the functional F(.) defined in (4.4), we have F(\(r) < F(0) = 0 
so that 

k 2 r- D 2 I |sinh(J(x,0 + *(x,0)'t<fa 

JD\K(t) 

< k 2 r~ 2 I cosh(^(x, t) + $(x,t))dx 

JD\K(t) 

.-2 



< k 2 r~ z / cosh(^(x, t))dx 

*D\KU) 



/_ 

Jd\ku) 

+ [ (Sy 0 fa(x,t)ds - [ (p$)(x,t)dx, 

J,)K(t) JD 
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and 



min(fc|.*2) / Wfa*,t)\ 2 dx 

JD 

,2 



1Ip(..OII6, 2 , d *iM..mj t2 , D IM-.'>HS,2.ajr(/) 
2*1 2 2€ 2 



+ g2H(yD^)(../)llg, 2 . ajf(f) 

" >\AT(0 



+ * 2 ro 2 ll*lloo.D / cosh($(x,OMx 



+ *3 + 2 

!!$(•, /)IIoo,dIIp(.,OIIq ! 2p €A\m,t)\\l 2D 

+ 2 U ' + 2 

By the trace theorem and Poincare's inequality there exist constants A.j = 
Xi(dK 0 ,D) > 0 and A 2 = k 2 (D) > 0 such that |fl>?(. , Ollo.2.aff(n 
< *iII?(.,0IIi.2,d and ||?(.,/)IIi.2,d < ^2l|V^(.,OBo,2,D. If we choose 
0 < e < min(*i,ifc 2 )/(2A. 2 A.f), t\ = <? 4 = Xfe, e 2 = 6 3 = 6 > 0 we have, 
for C\ = Ci(^i, fc 2 , Xi, X 2 ), C 2 = C 2 (/;i, k 2 , Aj, X 2 ), 



/, 



|VVr(X,OI^X 

D 

< C, (||p(. , /)|| 2 . 2<D + , Oll5.2.aiT(n) 0 + . OIIL.d) 
+ C 2 * 2 r- 2 ||$(.,r)lloo,i> f _cosh($(x,f)Mx. 
Finally, Poincare's inequality gives us 

Wl^(!0.T;H^D)) ^ C*(|D| 1/2 ||p|| L o C(OT:L o 0(D)) + ||a||o. 2 .3Ar 0 + ||*|| L c« ( o. r: / / .V2 (D)) ) 
(1 + ll*ll£.~(DT;WD») + C*||$|k-(DT;L-(Z»)llcosh$||^ (D 7 . :Ll(D)) , 

where 

C* = max ^XiCj /2 , XiC^kJ'VjJ 1 , C max fc 2 }^ . 

Then 

II VMI ^(0,7"; //^(D)) - 

where M = M (C*, *I>, ||p||L°wa«><z»), l[cr || 0 ,2.aAT(/))- 

The following theorem is central in order to establish an uniform bound 
forF (m) . 
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Theorem 4.1. Let us assume hypotheses (i)-(iv) then the solution \J/ {m) of 
(2 A) related with {K (m \t),p (m \o (m) ) satisfies, for each t e [0, T], 

max (\W\ m \. , f)|| 3/aA *« w . . OII 3/2 . 2 . D \^) * C ' 

where C does not depends on m. 

Proof. We only need to obtain the bounds for $ (m) = {^[ m \ fy m) ). Let us 
consider the problems 

V . (k(x, /)Vv (m) (x, 0) = 0, xe D, 
ufW) = i;j m) (x,r), x e 3/T (m, (/). 

u 2 m) (x,/) = 0, x€3D, (4,5) 

k 2 — 2 —(x, t) - 0 = - * (x< 0. x € 

on on 

and 

V • (k(x, /)V/ (w) (x, /)) - 6 (x, / (m) (x, t) + i/ w) (x, t) + $(x)) 
= p (m) (x,r), x e D, 



(4.6) 



/ 2 (m) (x,/) = 0, x g 3D, 
dn dn 

Using a variational formulation analogous to that in (4.3) we obtain solu- 
tions (weak) v (m \ f (m) e Hq(D) of the problems (4.5) and (4.6), respectively. 
We observe that f {m) + v (m) satisfies (2. 1) (in the weak sense), from the unique- 
ness of the variational solution for this problem, we have \// im) = f (m) + v lm) . 

The Theorem follows from the lemmas below. 

Lemma 4.2. Under hypotheses (i)-(iv) the solution v (m) e Hq(D) of (4.5) 
has the additional regularity 

v\ m) (. , f) 6 H y2 (K (m) (t)), u 2 m) (. , t) e H 3/2 (D\K^(t)) , 
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for each m>0andt e [0, T]. Furthermore 

max (\\v\ m) \\ LO c {O T . H y2 (K <m) (t))) , ll^r ) |lL^(OT;// ? ' 2 (D\jF^>(7)))) - C ' 

where C does not depends on m. 

Lemma 4.3. Under hypotheses (i>— (iv) the solution f (m) € Hq(D) of (4.6) 
belongs to C°(D) and, for each t e [0, 7"], 

supsup|/ (m) (x,r)| < +oo. 

m xeD 



Furthermore, ft m) (.,t) € H 3/2 (K im) (t)), / 2 (m) (.,0 € H 3/2 (D\K^(t)) and 
there exists C > 0 mc/i /ta 

max(||/, (m) || L ^ ( o,r;//v2 (A: (m. (0)) , ll/2 (m) ll L ~ (0 ,r:// 3 / 2 (D\^^))) - C ' 
w/i^re C does not depends on m. □ 



Proof of Lemma 4.2. Let us consider m > 0 and t e [0, 7"] fixed and define 
v\ m) = *, vf\ # 2 m) = k 2 vf\ then v = (v,, v 2 ) satisfies (in the weak sense) 

A# m) (x,/) =0, x € D, 

M2 uf } (x, f) = /), x e dK (m \t), 

t£ m) (x,0 = 0, xe3D, (4 - 7) 

3tJ < '" ) dv im) 

— 2-(x, f) - -^-(x, 0 = - (x, 0, x € a* (w) (0, 

where fi 2 = k 2 ~ l and i = fcf 1 . Following [3 1 ], we seek a solution v" = (T^i , V2 ) 
in the form 
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for£ (m) (.,0 € H l (dK lm) (0h <P {m) {.<t) 6 L\dK im) {t)) and x (m) (- < 0 
H\dD). Here 

(S (m V m) )(x,0 = f G(x-y)<p (m) (y,t)ds(y), 
(D<n (m) )(x, /) = f dG z X ~ y) ! im) (j,t)ds(y), (4 



(see [6] for details) 



' 0 






— o 






0 




_ x (w) . 



where 



and 



o , 2 (-/ +D r*)-„(-> /+ <>*) ,r»teP) 

(Dj m) C (m) )(x, /) = (p.u. D (w,) C (w) )(x,/), jr€8tf (m) (0 



(D^x^Xx, 0 = (p.v. D 0 x (m) )(x, r), jc € 8D. 
For each m € N and t € [0, 7*], the operators 

5 (m) : L 2 (3A- (m) (/)) H\dK (m \t)) 
Q/ + Dj m ^ : H\dK (m) {t)) H\dK (m \t)) 

y- X -l + D{™M : H\dK (m \t)) -+ H\dK {m \t)) 

Q/ + Dj"' r ^ : H l (dK im) (t)) L 2 {dK (m) {t)) 

(- X -l + DTJ : H\dK (m \t)) -> L 2 O^ (m) (0) 



/ + £>; 
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16 WEAK SOLUTIONS FOR THE ELECTROPHORETIC MOTION 

(jl + DSfj : H l (BD) ^ HHBD) 

are continuous [8] with operator norms depending only on K 0 , D. Observing 
that 

HVu(. , Ollo.2.8irc><o + IK n • Vv >(- . 0 ' n|lo,2.ajr<W(f> 

is an equivalent norm to ||.||i,2,dD ( see [32], Definition 1.9) and using the hypo- 
thesis |x - y| > d, Vx € dK (m) (t), Vy € dD, a direct calculation gives us that 
the following trace operators 

K jm) (0D (m). H \ {dD) ^ H \dK (m \t)), 

y£ m \t)D im) : H l (dK im) (t))^ H\dD) 
y< m \t)S (m) : L 2 {dK {m) (t)) -+ H\dD) 



(3D) -> L 2 (dK {m) (t)). 



are bounded and compact, with operator norms depending only on Kq, D and d. 
The compactness follows from an analogous argument as in [7] (Theorems 1 .6, 
1.7 and 1.10) if we observe that the kernels are continuous. 

Now it follows as in the paper by Torres and Welland [31] that J\. (mri exist 
and is a bounded operator on = H l (SK (m \t)) x L 2 {dK (m \t)) x H ] (dD) 
to X {m) . The operator norm ||J^ (mrl || depends only on d, K 0 , D. Then 



C (w) (.,0 
x (w) C,0 



<||A (mr, || P(.,/)|| 0 ,2,ar(")(r) 

X (m) 

< C\\JA^ ||(||o% , t)\\ 0XdK(mHt) + ||?(. , OII3/2.2.D), 



where X (m) is equipped with the product norm. 

As observed in [31], the operator norms of the operators in (4.9) and its ap- 
propriate inverses depend only on the Lipschitz character of the domain, so that 
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there exist L, = Li(3tf 0 ) > 0, L 2 = L 2 (dD, BK 0 ) > 0, such that, V/ e [0, T], 

< L, max , OliAlWw I* W C . Ollo.2,a*<»<,)) 

<L 2 max(||C (m) (.,OII,.2 ll^ <W) (.,Ollo,2 llx (m) (.,/)lh.2.aD). 

Hence, using hypothesis (iv), we have established that 

max (jiv\ m) \\ L oc (O T . H 3/2 (K(t))) , l|i4 m) lk~(OT;// 3 / 2 (Dvi^))) - c > 
where C does not depends on m. □ 

Proof of Lemma 4.3. Let us consider m > 0 and t e [0,T] fixed. The 
variational solution / (m) of (4.6) satisfies 



-k\ f (V^ (m) • V/ (m) )(x, t)dx -k 2 f (V^ (m) • V/ (m) )(x, t)dx 

JK^Ht) J D\K {m Ht) 

= / (p (m) ^ m) )(x,t)dx + k 2 r^ 2 f sinh(^ ) (x,r) + $(x))^ (,n) (x,/Wx. 

JD J D\K {m) (t) 

Then the Young's inequality gives us that 

^ f |V/<">(», t)\ 2 dx + ^ f |V/ (w) (x, t)\ 2 dx 

< - / (p (wi) ? (m, )(x, /)</x - * 2 r~ 2 /" sinh (£ (/w) (x, f) + $(x))£ (,n) (x, /)Jx. 

Using a similar calculation as in Lemma 4.1, the bound established there for 
\\$ {m) t , 0lh,2.D and hypothesis (iii), we have \\f m (. , OII1.2.D < C, where C 
does not depends on m and r. Observing that D is a C 2 -domain, standard elliptic 
estimates (see Chapter 14, Theorem 2.1 in [21]) show us that /<"*>(. , f) € C°(D). 
Let us define 

/? (B,) (0 = sup \f (m) (x, t)\ , then 3R(t) = sup /? (m) (r) < +00 . 

In effect, let us suppose that /?(/) = +00, then there exist a subsequence f (ntk) 
related with (n°**\ v (w * ) ), f (m *>, such that R imk) (t) -> +00 with k -+ +00. 
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Let us choose x^ e D such that /<"•*> (x< m *>) = R (mk) (t), then there exists 
8 {mk) > 0 such that 

|/ (m * ) (x (w * ) ,/)-/ (m * ) (x,/)l <2~ mk if xe B(x (mk \8 (mk) )nD. 

Recalling that f (mk) \dD = 0, if we take k oo we have, by the bound 
Wf (mk \.J)hxD < C, V) "> 0- In this case ||V/^(.,0l€AJ> -> +oo, 
which contradicts the bound \\f (ntk) (. , 0lli,2,J) < C. 

If we set = 6(x, / (m) + v im) + $) + /0 {m) we have /* (m) (. , 0 e L 2 (D), 
for all / € [0, T] and m > 0. Extending h {m) to be zero outside of D and setting 
g w = c * we have Ag (m) (. , /) = /i (,n) (. , t) a.e. in D and g (m) {. , f) e 
// 2 (D);if 



(m), 



D\^""»(l) 



we have 



^(..Ola^o € H\dK {m \t)), g™(.,t)\ aD e H l (BD) 



(see the proof of Theorem B in [20]), while ^-(. , t) — 
in 3A: (m) (r). Hence the solution / (m) of (4.6) can be written as (see [6]) 



(.,0 = 0 a.e. 



#1 l3ff(/) 

0 

L ft 0 "'^ J 



where 



4 w) J 



D (m) ^Sim) 
D (m) fl 2 S 

and f\ m) = k x f[ m \ T 2 m) = k 2 f± m) . The operators in tf™ were defined in (4.9) 
and JA im) was defined in Lemma 4.2. The uniform estimate in the // 3/2 -norm 
follows in a similar way as in Lemma 4.2 if we get uniform H 3 ^ 2 -bounds for 
g (m) (. , 0 and //'-bounds for g[ m \. , t)\ aK(t) , gf\. , t)\ BD . From the contin- 
uous imbedding H 2 (D) c H 3/2 (D) and by the boundedness of the operator 
F : L 2 (D) H 2 (D), where Fh (m \. ,t) = G* h {m) {. , t) = g (m) {. , t) (see the 
proof of Theorem 1 in [32]) we have the bounds 

IIS (W) (.,')ll3/2,2.D 

S ^(H VF=^*2 r D 2 Sinh + ^ + *)C . OII0.2.D + \\p (m) (. , OIIO.I.D). 
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where X depends only on D. A direct calculation and the above estimate give us 

llsi w) (. , OlliAMW ^ *V m) (. . Ollo.2.Dllri ,n) (. , OhxBD 
< XWow^kir? sinh (/<»> + v (m) + $)(. , /)|| 0 .2.D + l|p (w) (. . OII0.2.0) 

where A/ = k'(K 0 , k), k" = k"(D, k). Hence we need only obtain an uniform 
estimate for 

II j ' d\-k^hTM~d 2 sinh (f (m) + »™ + $ )(- ■ 0ll 0 ,2,x>\?55(o • 
We observe that 

sinh (f (m) + v (m) + $) = sinh + i/ m) ) cosh ($) 

+ cosh (/ (m) + u (,B) ) sinh ($). 

The result follows using the fact that *I> e L°°(D), recalling that by Lemma 4.3 
we have sup||/ (m) (. , f)llc°<D) < + 00 and from tne estimate 

m 

l|cosh(vf>)(.,f)|| 02 D ^ < Cexp (Cll^r^- , 0113/2.2.0x^7777) (4 lQ) 

<C, 

(see [6] and Lemma 4.2). □ 

Corollary 4.1. fj ||F (m) (. , t)\\l 1D dt < C, where C does not depends on m. 

As remarked in the preceding section the solution of (3.4) is constructed as 
a limit of a sequence of the appropriate approximation solutions and depends 
on the respective convergence of the force term F (w) to F in L 2 ((0, T) x D) 3 . 
Below we show in detailed manner how to do this. We begin with a technical 
lemma. 

Lemma 4.4. Let us assume hypothesis (i)-(v) and consider \jf (m) , \J/ the solu- 
tions 0/(2.1) related with (K (m) (t), p {m \ o (m) ) and (K(t), p, a), respectively. 
Then, \\\l/ {m) (. , t) - VK. , 0lh,2.D 0, with m -► +00, for each t e [0, T]. 
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Proof. Let us set r) m = \[r (m) - \f/. The variational formulation for \Jr, ^ (m) 
gives us 

- ki f (V^-VifcKx.Od* 

JK< m Ht) 

-k 2 f (V?™ • V^ w )(x, t)dx + /* 0* 

= I {p m r) m ){x,t)dx + k 2 r^ 2 I sinh(? (m) (x,/) + *(x,r))^(x,f)rfx 

J D\K (m Ht) 



>D\K (m Ht) 

and 



•/*0) 



- *2 / (V?. V^Xx,/)^ + f (o)w m )(x,t)ds 

JD\K(t) JdK(t) 

= / (pTi m )(x,t)dx + k 2 r~ 2 sinh(J(x,/) + *(x,/))r/ m (x,r)</x. 

Jd Jd\ku) 



If we define = tf (m) (0 n K(t), B {m) (t) = D\(tf(0 U K™{t)) we 

have, after subtracting the above expressions, 

- *i /" |Vifc(x, 0\ 2 dx - k 2 f |V^ m (x, t)\ 2 dx + f (ayorimHx, t)ds 
JA (m) U) JB< m >(t) JdK im *{t) 

- f (vYMm)ix, Ods = f (p {m) - p)(x, t)tf m {*, t)dx 
JdK(t) Jd 

+ kir ° 2 L^ t) 2cosh (t 1 + ^ + $ ) (x ' ° sinh (^f 1 ^)^^' 

-k 2 f |Vj(x,r)||V^(x,0|t/x 

JK< m >(t)\A im) U) 

+ ki [ \V$ im Hx,t)\\Vri m (x,t)\dx 
+ k 2 f \vfr m \x,t)\\Vri m (x,t)\dx 

- k 2 r~ 2 I sinh (f (x, 0 + $(x, r)Wx, t)dx 

J Ki m Ht)\Ai m Hn 

+ k 2 r~ 2 f sinh W (m) (x, t) + $(x, /))^ m (x, f)</x. 
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Using Young's and Holder's inequalities and positiveness of the second inte- 
gral on the right side, we have 

*l / |Vr; m (x,/)| 2 ^x + ^ f \Vrj m (x,t)\ 2 dx 
< I (a (m Vo^)(x, t)ds - f t)ds 

\JdK^(t) JdK(t) 

+ |D| 1/2 ||(p <w) - p)(. , OIIo.oc.dIIW- , OII0.2.D 

+ ki\K(t)\A {m) (t)\ ]/6 \\(S7f)(. , /)|| 0 ,3,jr(ol|V^(. , t)hxKu) 
+ k 2 \K*\t)\Ato\ t )\ l ' 6 wH . Ollo,3,jf(-)(o»V»; m (. , Ollo,2,jr<*>(,) 

+ k\\K <,m) (t)\A im) (t)\ l ^ (, \\V\j/^ m) (. , Ollo.3,A' (m, (/)ll^ 7 ^ W1 ^ ' ')Ho.2,/: <m, (/) 

+ fc 2 l^(/)\^ (m) (0| ,/6 ||V^ W >(, , f)||o,3,*(»)ll V^C , t)\\o.2,KU) 

+ k 2 r^ 2 \K^\t)\A im) (t)\ ]/4 \\ sinh {$ + *)(. , f )ll 0 ,2.*<«><i)ll»fo.(- , 0llo,4,*<"«>(0 
+ k 2 r- 2 \K(t)\A< m) (t)\ l/4 \\ sinh (£ ( "'> + *)(. , Ollo,2.« m (,) , 0 II 0.4, *„</)• 
Now, following [11], 

sup (\M im) (t) - M(t)\ + |M (m) (0 - 
f€(0,D (411) 

< C r \v™n™ - V p U\ L oo {0tT;LHD} y 0, 

with m -> +oo, since v ( p m) u {m) converges to v p u in C([0, T]; L 2 (D)) 3 . Hence 
\K im) (t)\A {m) (t)\, \K(t)\A {m) (t)\ 0, w -* +oo. We observe that, by the 
Theorem 4.1 and the Sobolev imbedding H i/2 (K {m \t)) c L\K im) (t)), 

l|V^ (m) (.,Ollo,3.^>(o < C|V* 1 w (.,0li A2 ,*w W 

< ll^, (W) (.,0ll3A2,^>(0 

< c, 

where C does not depend on m. Similarly, using the Sobolev embedding 
H l (D) c L 4 (D), Lemma 4.1 and hypothesis (iii) and (iv) we see that uniform 
bounds to \\ri m \\ 0 ,2,D, l|i?mllo.4,D are available. From hypothesis (v), ||(p (m) - 
P)(- . OIIo.oo.d -+ 0. Theorem 4.1 and similar argument as in (4.10) give us that 

|| sinh + $)(. , f)|oAa« w < Cexp (CH (m \- , 0 II 3/2.2, ) 

5 C, 
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where C does not depend on m. 
Using the rigidness of the particle we have 



/ {a (m) Yo r) m )(\, t)ds(\) - / (SyonrnHx, t)ds(x) 

JSKtoHt) JdK(t) 

= [ a lm) (M (m) (t)y, t)Yo(T}m(M {m) (t)y, t) - r) m {M{t)y), t)ds{y) (4.12) 

f Yon m (M(t)y, t){a (m) (M {m \t)y, t) - a(M(t)y), t)ds(y). 
JdKo 



Now, we introduce, for each m, the usual CJ°-regularization {n^} for n m such 
that Win™ - n m )(. , OIIi.i.d -> 0 when h -> 0. Observing that 



llyo^7^ ) (^/ (m) (Oy, 0 - ftqgkJfCtij, Ollo.2,a*„ -+ 0 with m -> +00 

the first integral in (4.12) tends to zero with m — ► +00, as can be seen using 
the usual trace theorem. Analogously we have the same result for the second 
integral in (4.12), recalling that a = ^f-a + (k\ ^ - *2i^). the hypothesis 
(iv) and regularity $(. , t) e // 3/2 (D). □ 

Theorem 4.2. Let us consider (p (m) e //'((O, T) x Df such that (p (m) -> <p 
strongly in C([0, 7*]; L 2 (D)) 3 and <p € #'(((), 7) x D) 3 . Then 

lim / / (¥ (m) • <p (w) )(x, r)£/xrfr = f [ (F • <p)(x, r)rfxrfr, 

for all t e (0, r). 

Proof. Recalling that = £- e (p ( 2 m) + r^ 2 * 2 sinh (^ m, ))(V^ m, )7 DxF= ^, 
we can write 

F (m) ■ (p m - ¥ • <p = Ci / DX ^(v> (m) " p) ■ V* w sinh (^>) 

+ C,<p • (/ DX ^7)V^sinh(^) 

- /^^V^^'sinh^^')) (4.13) 

+ C 2 «> • (l Dxm pV* ~ I D \J^/ m) V* (m) ) , 
where C, = and C 2 = 
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In what follows we shall estimate the integral of each term above. 
/ / 1 D\i^HT)^ ~ f) • sinh (* W »<* 

< \\<p (m) - <p\\ C( [0,T];LHD))) 

jf |V*«(. , r)ll 0 .3, Dx ^ll sinh <*«)(. , r)|| 0 6 . D ,^^r 

< C||^ (m) - p||(C[o.ru 2 (D))3 -* 0. 

with m +oo. Here we have used the Holder's inequality, the embeddings 
H l/2 (D) C L 3 (D)and// ,/2 (D) C f/ 3/2 (D), Theorem 4.1 jointly with a similar 
estimate as in (4.10) and the convergence \\(p [m) — <p||(cno.ri;L 2 (D))) 3 ~* 0. 

jC /d ( f ' ( W> v * Si " h W " V^7i V ^ <W) sinh (x ' t)</ju/t 

< f \[ \<p(x,r)-(ViJf sinh (^)-V^ (m) sinh (^ (m) ))(x, r)|</x 

/ 

+ / |?(x,t) • (V^ (m) (x.r))sinh(^ (m) )(x.r)Mx|^r. 

J/r(r)\Ar"")(T) J 

Writing 

Vt/r sinh (VO - V^ (w) sinh (\J/ {m) ) 
= sinh (V0(V^ - V^ (BI) ) + V^ (m) (sinh ftfr) - sinh 



(4.14) 



+ / \<p(x, T) (Vxf/(x, r)sinh(^(x. r))\dx 

Ht)\K(x) 



we have 



[ [ (p(x. r)(V^ sinh (xj/) -Vxl/ {m) sinh (^ (m) ))(x, r)<*x</r 
./0 JB {m Hr) 

< I || sinh (^r(. , r))|| 0 4 #<m)( r ) 

0.4.fl (m, (r) 

||<V* - V*<"»>)<., T)|| 0 2 .JM^)* (415) 

+ /' l*C • r )Ho,6,fl<«>(..r) • T )ll0.3.B<-( r) 

» 0 

||(sinh (+) - sinh , r)\\ 0 2tB ( mUx) dT. 
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Now, observing that if 

^(«) _ jff > o, cosh 2 (^ (m) + 0(xff - \J/ (m) )) < cosh 2 (t/r) and if 
Vr (m) - \/r < 0, cosh 2 (^ (m) + - xf/ (m) )) < cosh 2 (ifr - \ff {m) ) , 

we have, by Schwarz's inequality 
||(sinh (V) - sinh (V (m) ))(. , O\\ 2 0 2J}{mHt) 



<C\\W-xl, im) )(.j)\\ 2 OAJ) -+0, 

where we have used a similar estimate as in (4.10), Theorem 4.1 and Lemma 
4.4. Hence the dominated convergence theorem, Theorem 4. 1 and Lemma 4.4 
gives us that the terms in (4.15) tends to zero with m -► +oo. Passing the 
terms in (4.14) to the limit m -> +oo, using (4.1 1) and Theorem 4.1 we obtain 
the desired result for the second term in (4.13). The convergence of the others 
terms in (4. 1 3) to zero follows from a completely similar way, using hypothesis 
(iii), (iv) and (v), Theorem 4. 1 and Lemma 4.4. 'C 
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Abstract. We study the Riemann problem with forward combustion due to injection of air into 
a porous medium containing solid fuel. We neglect air compressibility and heat loss to the rock 
formation. Given initial reservoir and injection conditions, we prove that there is a unique time 
asymptotic wave sequence for combustion with complete oxygen consumption. The sequence 
consists of a region of unburned air at injection temperature, a warming discontinuity, a hot region 
with unburned air, a combustion wave and a region with burned air and unburned fuel at the 
initial reservoir temperature. The waves have very different speeds, and therefore they cannot 
be detected in laboratory experiments that focus on the combustion wave. However, they should 
occur in field scale. 

By introducing a cut-off in Arrhenius' law, the reaction rate vanishes at reservoir temperature, 
and two types of wave sequences are possible. One was already described. The other occurs for 
incomplete oxygen consumption. In this case, the wave sequence contains another wave, i.e., there 
is another region ahead of the combustion wave containing incompletely burned air at reservoir 
temperature, and a gas composition discontinuity that moves fast. However, instead of a unique 
solution for each Riemann data, there is a one parameter family of wave speeds . nd strengths. 
This multiplicity of solutions may to be due to the cut-off. 



#617/04. Received: 06/X/04. Accepted: 08/VII1/05. 



Copyrighted material 



28 WAVE SEQUENCES FOR SOLID FUEL 

Mathematical subject classification: 80A25, 76S05 

Key words: coflow filtration combustion, porous media, Riemann solution. 



1 Introduction 

Air injection and in-situ combustion have long been considered as a potential 
technique for displacement and recovery of heavy oil reserves [7, 22]. Opera- 
tional advantage of this thermal recovery technique is the abundance of injection 
gas independently of location. It utilizes heavy and immobile components of the 
crude oil as fuel producing in-place heat necessary for the recovery of upgraded 
crude oil. 

Despite the advantages and a long history, only a small fraction of the total 
thermal recovery utilizes this technique. Some reasons are technical, such as 
the possibility of front extinction and the necessity of (re-)ignition for sustained 
propagation within in-situ combustion in the presence of external heat losses [1]. 
Thus mathematical analysis of this problem is important to predict these events. 

A large number of studies on the structure of the combustion front have been 
reported since 1950s, see [1, 2, 3, 4, 5, 6, 8, 20, 21, 23, 26], for instance. 
These studies did not take into account other waves that occur in the combustion 
problem. As there is interaction between the combustion wave and other waves, 
this paper focuses on the solution of the Riemann problem with combustion 
taking into account all possible waves. 

We assume that upstream processes during the in-situ combustion have already 
generated a stationary homogeneously distributed fuel. Burning of this fuel is 
the subject of the paper. A bimolecular reaction is assumed to take place be- 
tween the injected oxygen and the solid fuel, hence the reaction region behaves 
as a source of heat as well as a sink for both of the reactants. We consider 
uniform flow, transport and reaction of injected air in porous medium of length 
/, with the reaction rate based on the Arrhenius' law. We neglect external heat 
losses. We develop simplified theoretical models for forward (co-flow) combus- 
tion under varying boundary conditions and obtain the wave sequences in the 
Riemann solution. The formulation of the problem is the same as in [4], while 
the nomenclature follows [1]. 
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The paper is subdivided as follows. In Section 2 we present the mathematical 
model. Theory on systems of conservation laws suggests that, besides the com- 
bustion front, there could be rarefaction, shock or contact waves. In Section 3 
we determine the characteristic speeds and prove that in this model there are 
neither rarefaction nor shock waves, only contact discontinuities. In fact, an 
examination of the Rankine-Hugoniot conditions shows that all noncombustion 
discontinuities are just the contact waves already found. In Section 4 we intro- 
duce the combustion front as a traveling wave of the evolution system derived in 
Section 2. In this section we also discuss the Rankine-Hugoniot jump conditions 
for the combustion wave. In Section 5 we prove that for the physical parameters 
considered in this paper, there are no resonances between waves, as opposed to 
the case analyzed by Aldushin et al. in [4]. As a consequence of the absence of 
resonances, we obtain two distinct temperature relationships for the combustion 
front, which we call the hot upstream and the hot downstream combustion cases. 
In Section 6 we determine the ranges of parameters in which these combustion 
waves may exist. In Section 7 we present our main results. For the hot upstream 
combustion we prove in Theorem 7.1 that the wave sequence of the Riemann 
solution is uniquely determined in case of complete consumption of oxygen. 
In Theorem 7.2 we prove that, for incomplete consumption of oxygen (using 
a modified version of the Arrhenius' law), there is another wave sequence as a 
solution. However, this solution is not unique, rather it is a one-parameter fam- 
ily of solutions. We also show that the hot downstream case does not occur in 
this model. In Appendix A we present tables with typical values of parameters, 
constants and nomenclature used through the paper. 

2 Formulation of the problem 

We assume that air is injected at the leftmost part of a porous rock cylinder 
containing solid fuel, so that all wave propagation is one dimensional. Balance 
equations are written for the total energy, the total gas mass, the oxygen mass and 
the fuel mass. For the latter, we define the fuel density per total volume p / and 
introduce the extent of conversion depth, n(x,t) = 1 - p f (x, t)/p° f (p° f is the 
initial fuel concentration), such that rj = 0 corresponds to complete availability 
of fuel (denoted by superscript o) and n = 1 to the complete lack of fuel (the 
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latter may occur because the fuel was never present or because it was completely 
consumed). The primary dependent variables are the temperature, f (x, t ), the 
oxygen concentration in terms of mass fraction, Y(x , t ) and the fuel conversion 
depth t)(xJ). The gas density p g (f , p) is expressed by the ideal gas equation 
of state in terms of temperature and total gas pressure p(xj). 

In formulating the conservation equations we make the following assump- 
tions: the pore space and the solid matrix are in thermal equilibrium so that a 
one-temperature model is used for the energy balance; heat transfer by radiation, 
energy source terms due to pressure increase, and work from surface and body 
forces are all negligible; the ideal gas law is the equation of state for the gas 
phase; thermodynamic and transport properties, such as conductivity, diffusiv- 
ity, heat capacity of the solid, heat of reaction, etc., all remain constant. We also 
neglect heat loss to the surrounding rock formation. We assume that gas heat 
capacity is negligible with respect to the rock heat capacity. The heat loss will 
be taken as zero as we study only the adiabatic case in this work. We assume 
that the pressure changes within any wave are negligible compared to the pres- 
sure drop across the system, so that in the ideal gas law and in other physical 
properties the pressure appears as a constant. Under these assumptions the di- 
mensional form (indicated by the superscript tilde) of the energy balance, the oxy- 
gen mass balance, the gaseous phase mass balance and the combustion kinetics 
equations are: 

d(c s Psf) B(c e p e vf) ~d 2 f 
(1 - 0) + = A— + Qp'W, (2.1) 

J(p,Y) f d(p g v?) d ( d?\ 

t-^r + -5- = Dm Vx (*im ) ~ {22) 

d(p g ) d(Pgv) . „ 
dr) 

Of = W, (2.4) 

where W is the reaction rate. The equations (2. 1 H2.3) were derived utilizing 
Darcy's law. The four equations (2.1)-(2.4) correspond to the three primary 
unknows T, Y, rj and the secondary unknow v, which is the volumetric flow 
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rate of the gas phase. In the above, c, denotes the average specific heat capacity 
of species / (gas or solid) at constant pressure, p, is the volumetric density of 
species i, X is the effective thermal conductivity (k s =X g ), Q is the heat released 
due to reaction (per unit mass of solid). Changes in the porosity 0 are considered 
to be negligible so that </> is constant. D M is an effective diffusion coefficient 
in the gas phase (D=D w /</>, where D is the molecular diffusion coefficient), 
while £l = yM a /Mf and ji gp = y gp M gp /Mf are mass- weighted stoichiometric 
coefficients for oxygen and for combustion gaseous products, respectively. M 0 , 
M f and M gp are the molecular weight of oxygen, fuel and gaseous products, 
respectively. The net gas mass production is determined from \x g = ji gp - /i, 
so that positive or negative sign for fi g correspond to net gaseous phase mass 
production or consumption, respectively. We will assume ix g > 0. For the rate 
of reaction, we use the second order law of mass action and the Arrhenius' Law: 



with activation energy E and pre-exponential factor k„. 

We use the ideal gas equation of state pM g = p g RT, where R is the universal 
gas constant, while M g and p g are the effective molecular weight and the effective 
density of the gas phase, respectively. We neglect the effect of changes of M g 
due to the reaction by approximating M g by a constant. 

Non-dimensionalized combustion equations. We introduce dimensionless 
space and time variables. To bring out the internal structure of the combus- 
tion wave, we introduce convenient variables and parameters that are defined in 
Appendix A. We scale the length by l*=a s /v' and the time using /*=/7u\ where 
v' is the injection velocity and <x s the effective thermal diffusivity. We introduce 
the scaled temperature 6 = f /7b, which means that the reservoir temperature 
corresponds to 0 O = 1; we also introduce the nondimensional gas density p, 
which is the gas density divided by the gas density at reservoir temperature. 
Thus the equations (2.1)-(2.4) are transformed into four dimensionless balance 
equations (2.6)-(2.9), while the dimensionless ideal gas law becomes Eq. (2. 10): 



W = k 0 e~ E/RT Y(\ -n), 



(2.5) 



d$ d(apv$) d 2 0 



(2.6) 



dt dx ' dx 1 
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"IT + IT = T e Vx{ p ^)~ (2J) 

dp d(pv) 

I = *• 

p9 = \, (2.10) 

where 

<t> = ae- y/G Y(\ -n). (2.11) 
The domain of the dependent variables is given by 

0>O, 0 < r < 1 , 0<f/<l, u>0. (2.12) 

From now on, we write 

The values of the parameters a, q, <f>, L e , n, fig, a and y are given in Table 1 
of Appendix A. 

3 Non-combustion waves for the model in hyperbolic framework 

In the absence of combustion, the source terms representing mass transfer or 
sensible heat generation containing the factor O vanish on the right hand sides 
of system (2.6M2.10). Of course <t> vanishes for Y = 0 or rj = 1. We focus 
in the waves for large times and long distances for which the second derivative 
terms are negligible. Let us first consider smooth solutions, so we can expand 
the derivatives in the remaining terms in Eqs. (2.6-2.9), manipulate Eqs. (2.7), 
(2.8), use Eq. (2.10) to eliminate p, obtaining: 

de dv 

U +a Y x = 0 ' (31 > 

8Y dY 

0— + u— =0, (3.2) 
dt dx 

(e \ do de rt 
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= 0. (3.4) 

dt 

The characteristic speeds of system (3.1M3.4) in increasing order and the 
corresponding characteristic vectors are: 

X"=0, (0,0,l,0) r , (3.5) 

0+a<f) \ 9+a<f)J 

A K = ^, (0,l,0,0) r . (3.7) 

It is easy to see that all characteristic speeds are constant along the integral 
curves defined by the corresponding characteristic vector fields, which means 
that all of them are associated to contact discontinuities, hence they satisfy the 
Rankine-Hugoniot conditions for (3. 1)— (3.4), [25]. The characteristic speed k* 1 
corresponds to an immobile discontinuity along which only n varies, X 6 corre- 
sponds to a thermal discontinuity along which 0 and v vary and k Y corresponds 
to a gas composition discontinuity along which only Y varies. 

4 The combustion wave 

We proceed next with the study of combustion wave propagation. The combus- 
tion wave encounters unburned states ahead and leaves burned states behind it. 

We look for combustion fronts as steady traveling waves of system (2.6H2. 1 1 ) 
with propagation speed V > 0 by setting x = x - Vt and t = t. In these moving 
nondimensional coordinates, after using Eq. (2. 10) to eliminate p, the equations 
(2.6)-(2.9) read, after writing x -> x : 



d 

d7x 



d_ 

Tx 
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In the adiabatic case studied here, reaction at the back of the combustion zone 
ceases due to complete lack of fuel, i.e., n b =l at -co. This is fuel-deficient 
combustion. (Thus, there can be no second combustion wave behind the first 
one.) Since gas is injected we have Y b = 1 . The temperature and the velocity 
behind the combustion need to be calculated. Thus the boundary conditions 
behind this combustion front are generically given by: 

9 = 9 b > 0, Y = Y b = 1, n = n b = 1, v = v b > 0; x -oo , (4.5) 

where the superscript b means burned. 

Ahead of the combustion zone we consider that there is abundant fuel, i.e., 
rj" = 0 at +oo, under two distinct conditions surrounding the combustion front. 
The first condition is the complete oxygen consumption case where we set K" =0 
at +oo. The second condition is the temperature-controlled case where we 
set 9" = 1 at +oo, and Y u have to be determined. The study of this second 
case is of special interest in petroleum applications, as it represents a situation 
where oxigen advances faster than the combustion wave, and the breakthrough 
of oxigen at the producing well causes safety hazards and other operational 
difficulties. These conditions can be summarized in (4.6) and (4.7), respectively: 

9 = 9 U > 0, Y = Y u = 0, n = rf = 0, v = v u > 0; x -+ +oo, (4.6) 

9 = 9 U = 1, Y = Y U > 0, n = rf = 0, v = v u > 0; x -> +oo . (4.7) 

where the superscript u means unburned. 

Integrating equations (4. 1)— (4.3) from x to +oo once, taking into account that 
n u = 0 at x = +oo and re-ordering, we get: 

Vi e-9 u ) + qn, (4.8) 



dY 



= L e 9 Q(v - <f>V)Y - l(u" - <t>V)Y u - nVn^j , (4.9) 



dx 

i {v _ 0V) _ l (v u _ 0V) + ^ Vr} = o, (4 . 10 ) 
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V-T = —$>■ (4.11) 
dx 

Next, we substitute the value of v given by equation (4.10) into (4.8) and (4.9), 
and substitute the value of <f> given by Eq. (2.11) into (4. 1 1 ) to obtain the reduced 
system: 

Y x = - ((£ - >)* + v ~ - v & -° u - «')■ < 412 > 

S = m ((f " + v ) y -^-tM- 0V)r )' (4 - 13) 

Sj* «y ( l_„ )e -r/». (4 .14) 
dx V 

We warn the reader that the temperature-dependent exponential on the right hand 
side of equation (4.14) will be modified when the temperature-controlled case 
given by the boundary condition (4.7) is studied. 

The existence of solutions of system (4. 12)— (4. 14) for the boundary conditions 
(4.5) and (4.6), or (4.7), connecting the burned state U b = (0*, Y b , n b \ v b ) at 
x -> -oo to the unburned state U* = (0 U . Y", rj"; v u ) and x -> -hoo is studied 
in a separate work, [9, 10]. Such solution represents the profile of a traveling 
wave connecting the burned state to the unburned state. For large times, this 
traveling wave is regarded as a thin combustion shock, [11, 15]. The Rankine- 
Hugoniot condition relating the burned and unburned states is studied in the next 
subsection. 

4. 1 The Rankine-Hugoniot equations for the combustion wave 

For the boundary conditions given in (4.5) at -oo and in (4.6), or (4.7), at +oo, 
taking into account that 86 /dx, dY/dx and dn/dx tend to zero as x tends to 
±oo, the R.H.S. of equations (4. 1 2)-(4. 1 3) become: 

° (t£ " ((!£ _ 0 0 + v " v ") " v(eh ■ ° u ~ q) = ° ' (4 - i5) 
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and Eq. (4.14) becomes a trivial identity. Taking x -» -oo in Eq. (4.10) we 
obtain a separate equation for v b : 

v h -<l>V= (J;(v''-4>V)-ii g V^e b . (4.17) 

The solutions of the Rankine-Hugoniot equations (4.15H4.17) together with 
the equation <J> = 0 determine all pairs of left and right states that are equilibria 
of system (4.12H4.14), i.e., that have the potential of representing burned and 
unburned states associated to a combustion wave. 



4.2 Geometric analysis of the Rankine-Hugoniot equations 

Since X Y is the particle speed of gas it is convenient to introduce the scaled 
combustion wave speed V M (V, v u ) = V/k Y (U u ) given by: 

<bV 

V u = — . (4.18) 
V 

Here the superscript u in V reminds the reader that V was scaled by v u . Equa- 
tions (4.15) and (4.16) give 0 b and Y" in terms of V, v" and 6", or in terms of 
V" and 0" as follows: 



e b = 



e b = 



«0 U +q+a<t>)V -av u )6 u 
((l+a^)6»"+a0) V -av u ' 

((6 U +q +a(f>)V u -a(f>)0 u 
((1 +ati g )d u +a<p) V -a(p' 



or 



(4.19) 



v" - (0 + (ji g + 11)0") V 
Y = — , or 

v u -4>V 



Y u = 



<P - (0 + (n g + nW) V 1 



(4.20) 



0(1 - V u ) 

For a fixed value of 9", Eq. (4.19b) represents 6 b in terms of V" as a hyperbola 
drawn schematically in Fig. 4. 1 , constructed for parameter values given in Ap- 
pendix A. The hyperbola has vertical and horizontal asymptotes respectively at: 



ad) h (0 U + q -I- a<b)9 u 

d {\+ati g )6 u +ad>' a (\+av g W u +ad>' 
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b 

e 




Figure 4.1 - 6 b as a function of V u given by (4.19b). 



The intersections of this hyperbola with the horizontal and vertical coordinate 
axes occur at: 

G h = 0\ and V" = — . (4.22) 

6 U +q+a4> 

Since afi g < q notice that Eqs. (4.21) and (4.22) lead to 0* > 0 U and V„ M < 
V l J < 1, respectively. 

Remark 4.1. The constants VJ and 9 b p = 0 b {V u p ) in Fig. 4.2 are related to 
Eq. (4.20b) and will be explained in the text. The constant 9 M is defined by 
Eq. (5.1) and V U M = V u (d M ) is obtained from (4.27). 

Now let us consider Eq. (4.20b) for a fixed value of 0" . This equation represents 
Y u in terms of V" as another hyperbola, which is drawn schematically in Fig. 4.2. 
This hyperbola has vertical asymptote V" = 1 and horizontal asymptote at 
Y " = <f> + 0 u ([jl + fA g )/<t>, which has no physical meaning; see the horizontal 
dashed line in Fig. 4.2. 

Comp. Appl. Math., Vol. 25, N. 1, 2006 



Copyrighted material 



38 



WAVE SEQUENCES FOR SOLID FUEL 




Figure 4.2 - Y u as a function of V" given by (4.20b). 



The intersections of this hyperbola with the horizontal and vertical coordinate 
axes occur at: 

V" = and r = 1 . (4.23) 

' 0 + (/i + ^)0" 

Knowing the value of V£ we define the scaled temperature value: 

0 b = 6 b (Vp, where the expression of 6 b is defined in Eq. (4.19b). (4.24) 
Notice that Vp = 1 / YJj < 1 . 

Remark 4.2. For the quantities in Fig. 4.1 and 4.2, corresponding to the pa- 
rameters in Table 1 , the inequalities 0 < V£ < V£ < < V p < 1 hold. 

We have expressed 0 b and Y u in terms of V". Now, let us also express 0 b in 
terms of the scaled combustion speed: 

V b = V/k Y (U b ) = <f>V/v b ; (4.25) 
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here the superscript b in V b reminds the reader that V was scaled by v b . 
Using Eq. (4.17) we obtain: 

v u = jl (e u v b + {{0 b - e u )<t> + fi g e u e b ) v^j . (4.26) 

Alternatively, instead of obtaining an expression for 9 b as a function of V" and 
9 U as in (4.19b), Eq. (4.15) could be utilized to find V" as function of 6 h and 6": 

a<b(9 u - B b ) 

V = _L__ (4 27) 

(0 tt ) 2 - {{9 b -q)- a(<f> - ii g e b ))9" - a<f>9 b ' 

Similarly V b can be written in terms of 9 b : 

, ad){9 b -9 u ) 

yb _ yv J m 28) 

" { e b ) 2 - (($" +q)- a(<f> + fi g 9"))0 b - a4>0" ' 




Figure 4.3 - V b (0 b ) defined by equation (4.28). 

Fig. 4.3 shows schematically V b (9 b ). The denominator in Eq. (4.28) is 
quadratic in 9 b and vanishes at the two real roots 9 b and 9 b , with 9\ < 0 and 
9 b > 9 U . Thus the graph of Eq. (4.28) possesses vertical asymptotes given by 

Comp. Appl. Math.. Vol. 25. N. 1, 2006 



Copyrighted material 



40 



WAVE SEQUENCES FOR SOLID FUEL 



9 b = 0f and 9 b = 9 b . The graph possesses also a horizontal asymptote given 
by V b = 0. 

We recall that we are interested in forward combustion, so we restrict our 
attention to the case V > 0. From Eqs. (4.18) and (4.25) it follows that V u > 0 
and V b > 0. Thus, only the portions in the first quadrant of the graphs in 
Figs. 4.1, 4.2 and 4.3 are considered. 

We conclude that the geometric analysis performed in this section yields the 
following restrictions relating burned and unburned states: 

o < v < v; , or vx < v < v; , (4.29) 
0 < e b < e u , or e b > e b > e u , (4.30) 

0<V*<1, (4.31) 

where V", Vj\ V n u , V;, 6 b and V b are defined by equations in (4.18), (4.21), 
(4.22), (4.23), (4.24) and (4.25), respectively. 

5 Wave resonances 

We can expect that changes of wave structure in the combustion Riemann pro- 
blem occur whenever two wave speeds coincide, either in the burned or in the 
unbumed regions. In this way we could have an interchange of the relative 
position of two waves. For the parameter values in Table 1 such changes are 
ruled out for physically interesting temperature range by the following theorem 
about absence of resonances: 

Theorem 5. 1. Let 6 M be defined by 

e M = —. (5.1) 

Then in the physical domain defined by Eq. (2.12): 

(a) the non-combustion waves have distinct speeds; 

(b) the combustion, the immobile and the gas composition waves have distinct 
speeds; 
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(c) consider a combustion wave following a thermal wave, i.e., V < X d {U"). 
In the limit as they tend to have equal speed, the temperature at the left of 
the superposed waves is 6 h = 6 M . 

(d) consider a thermal wave following a combustion wave, i.e., X e (U h ) < V. 
In the limit as they tend to have equal speed, the temperature at the right 
of the superposed waves is 6" =6 M - 

Proof. 

(a) From equations (3.5)-(3.7) we obtain explicitly that 

X n < X 6 < X Y everywhere, (5.2) 

which means that there are no speed equalities among non-combustion 
waves. 

(b) Since we are considering forward combustion there is no resonance bet- 
ween the combustion and the immobile wave. 

From Eqs. (4.29) and (4.31), we have that 0 < V" < 1 and 0 < V b < L, 
respectively. Since V" = V/X Y (U U ) and V h = V/\ Y (U b ), it follows 
respectively that 

V<X Y (U U ) and V < X Y \U b ) . (5.3) 

(c) From Eqs. (3.6), (4. 1 8) and (4.27), the combustion wave following a ther- 
mal wave tend to have the same speed if, and only if, V u = a <f>/ (9 U + a<f>), 
or 

(0" -0 h ) _ 1 

( 0 U) 2 _ ( (0b _ q) _ ai( p _ H e*>))9» - ~ e-+ad>' or (5 4) 

9 u {afi g 0 b -q) = 0 or 0 b = 0 M . 

(d) The proof is analogous to that of (c), using Eqs. (4.25), (4.28) instead of 
Eqs. (4.18), (4.27). This conclude the proof of Theorem 5.1. 
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Direct calculations with the value of the parameters in Table 1 yield 0 M % 
24.21264641. This corresponds to an extremely large temperature value above 
6000 K that is beyond physical interest in our case; even at temperatures lower 
than 6000 K we may expect that heat losses to the surrounding reservoir can- 
not be neglected. This value 6 M defines a maximum temperature value, below 
which there are no speed coincidences. Thus we focus our attention to the tem- 
perature range 0 < 6 0a/- 

So one conclusion of this section is that for physical values of temperature 
there are no resonances between waves in this model. 

Remark 5.2. The case where the thermal wave speed coincides with the com- 
bustion wave speed yields high temperatures and therefore maximal combustion 
efficiency. It is called superadiabatic combustion, which is desirable in other 
practical applications. The structure of such a combustion wave was analyzed 
in [4]. 

Since V" and V b are both positive, by analyzing the first denominator in 
Eq. (5.4) one can see that the sign of {G b - 0") determines if V < X 9 or V > X 9 . 
Thus, our final conclusion in this section is that the following two temperature 
relationships for the combustion front are possible: 

(A) Hot upstream (left) 

0 b >9 U > 0, for V > X 9 (U b ) and V > X 9 (U") ; (5.5) 

(B) Hot downstream (right) 

0 < G b < 6\ for V < X 9 (U b ) and V < X 9 (U U ) . (5.6) 

6 The admissible Rankine-Hugoniot locus for the combustion wave 

Cases (A) and (B) in (5.5) and (5.6) yield distinct possibilities for the ranges 
of 6 b , V b , V" and Y u corresponding to distinct portions of the branches in 
Figs. 4. 1 , 4.2 and 4.3. Such ranges are called the admissible states for combustion 
waves. In the next two subsections we will obtain the admissible ranges, which 
are fully characterized by the value of V": the first one is VjJ < V" < Vjj and 
the second one is 0 < V u < V". 
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6. 1 Hot upstream 

Since 9 b > 9 U (see Fig. 4.1) it follows that V u > Vf. As 0 < Y b < 1 (see 
Fig. 4.2) it follows that V u < V". 
Let us define 

V£ = V B (0 w ), (6.1) 

represented in Fig. 4. 1 and obtained from Eq. (4.27), where 9m is the maximum 
value for the nondimensionalized temperature 9 defined in (5.1). Let us also 
define 

V M = V b (9 M ), and V* = V*(#, (6.2) 

represented in Fig. 4.3, where the expression V h is defined in Eq. (4.28). Finally, 
let Y U M in Fig. 4.2 be defined as 

Y U M = Y U (V M ), where the expression for Y u is defined in Eq. (4.20b). (6.3) 

By comparing the graphs in Figs. 4.1, 4.2 and 4.3, and taking into account the 
inequalities obtained through the previous sections we arrive at the following 
Lemma, which characterizes the admissible Hugoniot locus for hot upstream 
combustion waves: 

Lemma 6.1. Let 9" be a fixed value of the scaled temperature in the range 
0 < 9" < 9 M . Then the inequalities 

V" <V U < V" , 
M - p (6.4) 

(see Figs. 4.1, 4.2 and definitions in (6.1), (4.18), (4.23)), 

hold if, and only if all of the following inequalities hold: 

9 b p < 9 h < 9 M , (see Figs. 4.1, 4.3 and definitions in (4.24), (5.1)); (6.5) 
0 < Y u < Y U M , (see Fig. 4.2 and definition in (6.3)); (6.6) 
V M <V b < Vp , (see Fig. 4.3 and definitions in (6.2) and (4.25)). (6.7) 

We notice that condition (6.5) implies that 9 b > 9 U (see Figs. 4.1 and 4.3). 
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6.2 Hot downstream 

Since 0 < 0 b < 6 U we see that 0 < V u < V„", where V n " was defined in (4.22). 
See Fig. 4.1. Let Y% in Fig. 4.2 be defined as 

y" = Y U (V»), (6.8) 

where the function Y u is defined by Eq. (4.20b). 

As 0 < V" < and 0 < 6 b < 0 U it follows that y„" < Y u < 1 and 
0 < V b < 1. See Figs. 4.2 and 4.3, respectively. 

As in the previous case, by comparing the graphs in Figs. 4. 1, 4.2 and 4.3, and 
taking into account the inequalities obtained throughout the previous sections 
we arrive at the following Lemma, which characterizes the admissible Hugoniot 
locus for hot downstream combustion waves: 

Lemma 6.2. Let 9" be a fixed value of the scaled temperature in the range 
0 < 6" < 9 M - Then the inequalities 

o < v it < v; , 

(see Figs. 4.1, 4.2 and definitions in (4.18) and (4.22)), 
hold if, and only if all of the following inequalities hold: 

0 < B b < 0", (see Figs. 4.1, 4.3), 
Y" < Y u < 1 , (see Fig. 4.2 and definition in (6.8)), 
0 < V b < 1, (see Fig. 4.3 and definition in (4.25)). 

7 Wave sequences in the Riemann solutions 

We have completed in Section 5 the proof that there is no wave speed coinci- 
dence in our combustion problem for temperatures of physical interest. We have 
also determined in Section 6 the admissible ranges of the left and right states 
of the combustion wave along the Hugoniot locus for the hot upstream and the 
hot downstream combustion cases given in (6.4)-(6.7) and (6.9M6.12), respec- 
tively. We will show in Section 7.2 that the hot downstream combustion case 
is unacceptable for this model. In the hot upstream combustion case we will 
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describe the wave sequence in the Riemann solution under the conditions (4.6) 
and (4.7) surrounding the combustion front, where the waves satisfy an extension 
of Lax speed inequalities. 

7.1 Hot upstream combustion case 

In the hot upstream combustion case the thermal wave precedes the combustion 
wave (X e (U b ) < V), there is no temperature change ahead the combustion 
wave, which means that 6" = 1 (reservoir temperature) and 6 b > 6 U = 1. Thus, 
generically the wave sequence in the Riemann solution consists of a (perhaps 
trivial) immobile fuel shock, a thermal shock with speed X 9 , a combustion front 
with speed V and a gas composition wave with speed X Y . We denote this se- 
quence of waves by means of the following convention: 

Ui zjt 2l+ JL> u» v o . (71) 

The state V = (0 1 , 1,0, v' ) denotes the injection conditions, U 1 = (0\ 1, 1, v') 
denotes an intermediate state in the burned region, while U b = (6 b , 1, \ ,v b ) 
and U u = (1, Y u , 0, v li ) are the burned and the unburned states surrounding 
the combustion front and U° = (1, 0, 0, i>°) denotes the reservoir conditions 
at production. The values of 6' and v' are given as boundary conditions, the 
values of 6 U and Y" are given by condition (4.6) or (4.7), but the speeds X°, V, 
X Y and the values of 9 b , v b , v" and v° have to be determined. 

7. /. J The hot upstream combustion in the oxygen deficient case 

As discussed above, we have that 9 U = 1 and according to the condition (4.6), 
which characterizes the oxygen deficient case, Y" = 0. We have the following 
theorem, which provides formulae for all the states as well as speeds for com- 
bustion and noncombustion waves in the wave sequence (7.1) for the Riemann 
solution in this case: 

Theorem 7.1. Assume that in the wave sequence for the Riemann solution 
of (2.6)— (2.1 1) there is a hot upstream combustion wave with left and right 
states satisfying (4.5)-(4.6) and d b > Q u = 1. Gi ven the injection conditions 
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U' = (0', I, 0, v'), with 6' > 0 and v' > 0, then the constant states and the 
speeds of all waves in the wave sequence (7.1) are uniquely determined. 

Proof. Inspecting Fig. 7. 1 we see that the speeds and the intermediate states in 
the wave sequence (7.1) are determined as follows. 

Substituting Y u = 0 in Eq. (4.20b) it follows that V u = V u p , where V u p is given 
in (4.23) with 9 U = 1 (see Figs. 4.1 and 4.2). From Eq. (3.6) and Eq. (4.19b) 
with 0 U = 1 and V u = V;, it follows that: 

X° = -^ and g» s l±iZffr±M 
0' +a(f> 1 - an 

Taking into account that the thermal wave is a contact discontinuity and Eq. (4.28) 
with 0" = 1 it follows that: 

v = — v and V = 



6" +a0 (6 b ) 2 - (1 + q - a{<t> + fi g ))9 b - a(f>' 

From Eq. (4.25) we obtain the combustion speed V; then we obtain v u from 
Eq. (4.18): 

v b V b „ 0V 

Since Y u = 0, the strength of the gas composition wave is zero, which explains 
why this wave is discarded in the wave sequence (7.1) as represented in Fig 7.1. 
Thus the proof of Theorem 7.1 is complete. 

7.1.2 The hot upstream combustion in the temperature-controlled case 

According to condition (4.7) we have 6 U = 1, but Y u have to be obtained. 

First of all we notice that the condition (4.7) cannot be analyzed as (4.6) since 
the unburned state U u is not an equilibrium of system (4.12)-(4.14), because 
<P does not vanish. However the exponential factor in (2.5) is extremely small 
at prevailing temperatures, so we use the following modified version of the 
Arrhenius' law, [27]: 

W(f) = k 0 e- Eb/mf - fo)) Y(\ -n), for f > f 0 , and 

1 . (7.2) 
W(T) = 0, forT<T 0 , 
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Figure 7. 1 - Regions separated by immobile, thermal wave and combustion waves in 
the Riemann solution. Values of 6, Y, tj and v in each region. 



where E b = E(f b - f 0 )/f b . 

This modification on the denominator in the exponential factor ensures that 
the reaction rate is equal to zero below the reservoir temperature To, which is a 
good approximation of the original formula in the Arrhenius' law. The value of 
E b is chosen so that at the burned state above reservoir temperature the exponent 
assumes the same values as in the original Arrhenius' formula. Thus in the scaled 
temperature variable the function O in (2. 1 1 ) is replaced by the approximation 

<& = cte~ r W b) Y(l -17), for 6» > 1 and <D = 0 for 9<\. 

With the above modified version of Arrhenius' law, the unburned state U u 
becomes an equilibrium of system (4.12)-(4.14) and the existence of traveling 
wave solutions can be obtained, [9]. 

Now we notice that since 6" = 1 the range for V" in Eq. (6.4) can be restricted 
further to 

V&<V U <V?, where Vf = — !■ (7.3) 

0 + /X + fl g 

is obtained from Eq. (4.23) with 0" = I. Consequently we get 

0? < e b < o M , 0 < r < y u m , v b M < v b < v b , 

where 0 b = e b (V l u ) and V b = V b (0 b ) are obtained from Eqs. (4. 19a) and (4.28), 
respectively. 
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We can prove the following theorem in this case: 

Theorem 7.2. Assume that in the wave sequence for the Riemann solution of 
(2.6)-(2. 1 1) there is a hot upstream combustion wave governed by version (7.2) 
ofArrhenius ' law with left and right states satisfying conditions (4.5) and (4.7), 
and B > 9" = 1. Given the injection conditions U l = (9' , I, 0, v l ), with 
9 > 0, V > 0 and any positive value of the scaled combustion wave speed V" 
in the interval [V^, V"] given in (7.3), then all the constant states and speeds 
of all waves in the wave sequence (7.1) is uniquely determined. 

Proof. Inspecting Fig. 7.2, we see that if the injection rate v' , the temperature 
9 l and the value of V" are given, then the speeds A. 6 *, V and k Y as well as the 
values of 9 b , v b , Y u , v u and v° are determined as follows. 
From Eq. (3.6) and Eq. (4.19) with 9 U = 1, we have: 

9 av i Qb (1 + q + a(t>)V u - a<t> 



r = and 9° = 



9'+a<p ' ' (1 +a^ 5 +<*</>) V u -acp ' 

Taking into account that the thermal wave is a contact discontinuity and 
Eq. (4.28) with 9 U = 1 it follows that: 

" {0 b ) 2 - ({\ + q) - a® + H g ))0 h - a* 

From Eq. (4.25) we obtain the combustion speed V; then we obtain v" from 
Eq. (4. 1 8). From Eq. (4.20b) with 9 U = 1 , we obtain Y u and from Eq. (3.7) with 
v = v u , we obtain the value of X Y : 

4> ' ~ V ' <hi-v) ' 0 ' 

Finally, since the gas composition wave is a contact discontinuity it follows 
that v° = v" , which provides a consistent solution of the system of conservation 
laws in the hyperbolic framework under the conditions (4.5H4.7) for left and 
right states to the combustion wave. This completes the proof of Theorem 7.2. 
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V 



Figure 7.2 - Regions separated by immobile, thermal, combustion and gas composition 
waves in the Riemann solution. Values of 0, Y, t? and v in each region. The latter region 
is actually very thin. 

Remark 7.3. We notice that, due to the fact that the combustion and the ther- 
mal waves are very slow, it takes a long time for these waves to separate from 
each other, while the gas composition wave (the extremely fast wave) separates 
from the others immediately. This should explain why such phenomena have 
not been detected in laboratory experiments, where transient rather than asymp- 
totic behavior is detected. 

Remark 7.4. As we have seen in Theorem 7.2, in this case the injection con- 
ditions together with the initial data are not sufficient to determine the Riemann 
solution uniquely, but they determine a family of Riemann solutions depending 
on the parameter V". This strange multiplicity of solutions for the Riemann 
problem may to be related to the modification of the exponential factor (7.2) in 
Arrhenius' law. 

Remark 7.5. There are two types of Riemann solutions for any given data 
if we use the modified version of Arrhenius' law: the regular solution for 
oxygen-controlled combustion and the one-parameter family of solutions for 
temperature-controlled combustion. 
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7.2 Hot downstream combustion 

Recall that in the hot downstream combustion case we have to restrict the pa- 
rameters to the ranges defined by Eqs. (6.9M6.12). Since 0 < y„" < Y u < 1 
the complete oxygen consumption case defined by Eq. (4.6) is impossible. 

Remark 7.6. We conclude that this model does not support case (b) of Aldushin 
et al. [4] that they call the reaction trailing structure. 

In the temperature-controlled case defined by (4.7), under the standard Arrhe- 
nius' law (2.5) the state ahead the combustion front is not an equilibrium of the 
system of ordinary differential equations (4. 12)-(4. 14). This fact can be verified 
by a procedure analogous to that employed in Subsection 7. 1.2. Thus the travel- 
ing wave technique is not applicable. On the other hand, if we use version (7.2) 
of Arrhenius' law we must set 9" = 1 and therefore we should have 9 b < 1 in 
the downstream combustion case. Thus drj/dx is alway zero, i.e., there is no 
fuel consumption near the unburned state U b and there is no connection between 
the equilibria U b and U*. 

8 Discussion 

We have seen that the original formulation of the Arrhenius' law allows us 
to describe mathematically the combustion front as a traveling wave for hot 
upstream combustion with complete oxygen consumption ahead the combustion 
front. In this situation we have proved that the wave sequence in the Riemann 
solution is uniquely determined by the injection conditions together with the 
initial data. 

For incomplete oxygen consumption, or temperature-controlled case, the 
original formulation is inappropriate to characterize the combustion front as 
a traveling wave, since the unburned state is not an equilibrium point of the 
associated ordinary differential equation system. In such case, in order to apply 
the traveling wave theory a change of the Arrhenius formula can be used so that 
the unburned state becomes an equilibrium. As a consequence of this change 
of the Arrhenius' law, we lose uniqueness of Riemann solutions. Although the 
wave sequence is similar to that for the complete oxygen consumption case, the 
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intermediate constant states and the wave speeds form a family depending on 
one parameter. 

It is hard to establish the existence of traveling waves, which is assumed in 
this work. Very often we are led to non hyperbolic equilibria of the associated 
system of ordinary differential equations in high dimensional spaces. Even 
in cases where the analysis can be done by the invariant manifold theory, we 
establish the existence, but without explicit formulas, which would be very useful 
[13, 17, 19]. Of course, an alternative approach to be used together the inva- 
riant manifold theory is the singular perturbation theory, also called "matched 
asymptotic expansion", in which explicit approximations of the traveling wave 
solution can be obtained, [9, 12 14, 16, 24]. This is the subject of a separate 
work, [10]. 
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Appendix A. Topical values, nomenclature and constants 

In Eqs. (2.6)-(2.12), we introduced the following variables and parameters 



x = 



v 




P-Po 




V 



(A.l) 



P = 



Pinj - Po 



(1 = 



p P * = 



Ppzp'i 

P>Y' ' 



Ps = 



~ n 

PsPf 



P'Y'' 



(A.2) 



a = 



(1 -0)c,p/ 
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q = 



Qp° f 

{\-<b)c s p s f 0 



K(Pinj - Po) 



K = 



h = 



ht* 

(\-4>)c s PsfoH 



(A3) 



<*s = 



(1 -0)C,p/ 



L e = 



a. 



y = 



RT 0 



a = koY'poU, (A A) 



where p Q corresponds to the initial gas pressure and is typically much larger than 
the pressure drop across the system. 



Physical quantity 


Symbol 


Value 


Total heat content of the porous medium 


<i 


1.0121 


dimensionless stoichiometric coefficients for oxygen 


1* 


205.8 


dimensionless stoichiometric coefficients for gaseous products 




68.19 


Lewis number (ratio of thermal and molecular diffusion) 


L e 


0.214 


Arrhenius number (dimensionless activation energy) 


Y 


23.69 


dimensionless reaction coefficient 


a 


0.027 


volumetric heat capacity ratio of the filtrating gas and matrix 


a 


6.13E-4 


porosity of the medium 


<t> 


0.3 



Table 1 - Typical values of dimensionless parameters. Source: [1], [2]. 
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Physical quantity 


Symbol 
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To 
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n 

Of 
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Abstract. Acid injection in porous medium is a process widely used for stimulation of 
petroleum wells and leads to the formation of highly conductive channels called wormholes. Two 
different transport-reaction models have been developed in Part I to describe the phenomenon at 
the core-scale. The possible existence of core-scale effective properties which appear in these 
models is discussed here on the basis of Darcy-scale numerical experiments. The advantages and 
drawbacks of one-equation and two-equation models are investigated by reference to averaged 
fields computed from Darcy-scale simulations. 

Mathematical subject classification: Primary 74Q20; Secondary 80A32. 

Key words: core-scale dissolution, acidification, effective properties, porous media, wormhole. 

1 Introduction 

The unstable dissolution of a porous medium leads to complex patterns which are 
difficult to model quantitatively [4]. Indeed, this dissolution process is coupled 
with the fluid momentum equation in an unstable way: flow velocity is higher in 
the largest pores, which generally produces faster dissolution processes. These 
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processes increase locally the pore diameter and this may in turn facilitate the 
acid transport to these large pores. These physical mechanisms can lead to the 
formation of highly conductive flow channels called wormholes. A model has 
been developed at the Darcy-scale by Golfier et al. [1] involving a local non- 
equilibrium dissolution equation. Numerical simulations have been performed 
for 2D and 3D configurations for both homogeneous and heterogeneous sys- 
tems and allowed to capture all the observed features in terms of dissolution 
regimes and optimum flow rate. However, a direct application of laboratory re- 
sults to the field scale is not straightforward since a direct Darcy-scale description 
would require a very fine grid. Therefore, a large-scale model is necessary. 



Darcy-scale Core-scale Field scale 




Figure 1 - The different scales of the problem. 

A first attempt has been made in Part I [2] at deriving a core-scale dissolu- 
tion model based on cross-sectional averages from a volume averaging method 
[5, 10]. We remind below the notations used at the core-scale and briefly sum- 
marize the different approaches developed in Part I [2]. It must be emphasized 
here that this study does not concern only stimulation methods by acidification 
of petroleum wells but the dissolution mechanisms in porous media in a more 
general way. As a consequence, we consider only the interaction of the water- 
acid mixture with the rock (single-phase flow) and the oil phase is not taken 
into account. Core-scale averaged quantities are defined in the core-scale av- 
eraging volume illustrated in Fig. 1. For example, the regional intrinsic acid 
concentration and the superficial regional velocity associated to the m— region, 
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respectively C* Am and V^, are defined as follows. 

C' Am = {C A Al = — ( C Ap dV (1) 

v uj JV m 

V; = {V # L = ^/^V,rfV (2) 

Flow equation. For the flow description at this scale, a classical Darcy's law 
has been obtained by averaging the Darcy-Brinkman formulation. 

v; = -— K*.(vp;- P/Jg ) O) 

/■*/) 

v • v; = o (4) 

where V£ and P£ represent respectively the core-scale averaged velocity and 
pressure and K* the core-scale permeability tensor. It must be pointed out that 
this average velocity is linked to the two regional average velocities and V*, 
with the following relation: 

v; = v; + v; (5) 

With regard to the transport and dissolution part, it was not obvious to apply 
an upscaling method which would lead to some core-scale equations valid in 
a general way. In fact, the value of the mass transfer coefficient appearing in 
the Darcy-scale model can strongly modify the acid transport behavior and the 
corresponding Darcy-scale dissolution pattern. Several approaches have been 
developed depending on whether the local mass equilibrium condition is verified 
or not. 

Local mass equilibrium dissolution: Two-equation model. If the value of 
the mass transfer coefficient is very large, we obtain two different regions at the 
Darcy-scale: a fluid region in the dissolved areas, and a porous region where 
dissolution has not occurred. This looks like a double-porosity system, and 
I this suggests the introduction of a two-equation model for which the wormholes 

i 
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(m— region) and the remaining porous matrix (//-region) are treated separately. 
The transport equations for the or— region are written as 

+ v - * vc ^ = v " ( D - vc ^) - a * c *<» (6) 

d<Pm _ ft * * 



a* a 



'a 



where & represents the stoichiometric coefficient of the reaction, p„ the rock 
density, a* the core-scale mass transfer coefficient and <f> m the wormhole volume 
fraction whereas the average concentration associated to the r/— region, CJ , is 
equal to 0. 

Local mass non-equilibrium dissolution: One-equation model. In the more 
general case of local non-equilibrium conditions at the Darcy-scale, the more 
complex form of the local equations did not allow to directly infer the core-scale 
transport equation. Nevertheless, it was possible to propose a general form for 
the macroscopic transport-reaction equations within a one-equation model [6, 
9], defined as follows: 

vC A ft 

+ V; • VCJ, = V • (DJ'.VC^) - a'CJ, (8) 
= -«*C% (9) 



Bt 0 



where e*p represents the core-scale porosity and the average concentration 
weighted by the porosity. 

While we have now decided on the form of the core-scale equations to be 
used in our model, a fundamental question remains: how to determine the differ- 
ent effective coefficients that appear in these core-scale equations? They could 
be calculated in principle by solving so-called closure problems. However, this 
requires the knowledge of the developed wormhole geometry, which is the result 
of the flow process. As a consequence of the fluid-porous medium interface 
evolution during the dissolution process, the relationship between the macro- 
scopic quantities may be historical. This problem exists also for the Darcy-scale 
model, and we have adopted in this case the assumption classically made in geo- 
chemistry of a direct relationships linking the different macroscopic properties. 
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For instance, permeability is a function of porosity, or the mass transfer coeffi- 
cient depends on the cell Peclet number and porosity. Is it possible to use some 
similar relations at the core-scale? 

In order to test such a possibility, rather than to solve numerically the closure 
problems with some non-representative boundary conditions, the dissolution pat- 
terns obtained at the Darcy-scale are used in this work to determine the effective 
coefficients. A series of Darcy-scale simulation is performed on a core of 25 
cm length and 5 cm wide and used to obtain, by spatial integration, some core- 
scale effective parameters. We remind briefly in a first part the various results 
extracted from Darcy-scale simulations and presented in Golfier et al. [ 1 ]. Then, 
the definitions associated with the introduction of the core-scale parameters are 
presented and the possible independence of these coefficients with respect to 
the dissolution history is discussed based on the numerical experiments. 

2 Darcy-scale numerical simulations 

Different Darcy-scale simulations have been performed on a two-dimensional 
domain by varying the injection flow rate Q, the initial acid concentration C 0 and 
the mass transfer coefficient value a. Two principal results have been extracted 
from the model. First, the simulations allow to capture the different dissolution 
regimes, those corresponding to local equilibrium as well as those with non- 
equilibrium. We can see in Fig. 2 from Golfier et al. [1] an example of the dif- 
ferent dissolution figures obtained numerically for different acid injection rates. 
All these figures represent the porosity field. The obtained dissolution patterns 
are remarkably equivalent to the experimental dissolution patterns. These results 
were used to draw diagrams of the transitions between the different regimes. 

Secondly, the performed simulations have confirmed the existence of an op- 
timum injection rate. This numerical optimum injection rate has been fitted 
with the experimental value through the proper choice of a single parameter, A, 
defining the mass transfer coefficient value through 

a (Pe ce u, £ P ) = Aa 0 (Pe ceU , Bp) (10) 

where or 0 is the correlation obtained from the "closure problem", presented in 
Golfier et al. [1]. This latter correlation expresses the mass transfer coefficient 
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(a) (b) (c) (d) (c) 




Figure 2 - Porosity fields representative of the dissolution structure obtained numerically: 
Example of dissolution patterns: (a) face dissolution, Pe = 8.32 x 10~ 4 , Da = 120;(b) 
conical wormhole, Pe = 4. 14 x 10 -3 % , Da = 24; (c) dominant wormhole, Pe = 1 .66, 
Da = 6.01 x 10 -2 ; (d) ramified wormhole, Pe = 83.2, Da = 1.2 x 10 -3 ; (e) uniform 
dissolution, Pe = 832, Da = 1.2 x 10~ 4 (Golfier et al., J. Fluid Mech., 457, 2002). 

as a function of the Peclet number, defined at the Darcy-scale Pe ce u, and the 
porosity ep. 

The calculations show that the optimum injection rate is linked with the Dam- 
kohler number value, called Da 1 . Also, this model allowed to study the effect of 
many external parameters, such as the variation of temperature or concentration, 
or the use of other reacting fluids, like acids in emulsion, which give lower 
diffusion effects. 

3 Macroscopic parameters and core-scale relations 

In this section, we discuss the determination of correlations between the dif- 
ferent core-scale effective coefficients from numerical experiments performed 
at the Darcy-scale. The definitions associated with these core-scale parameters 
are discussed below. We consider the domain represented in Fig. 1 . Averaged 
properties, that could be used in a ID core-scale representation are defined as 
follows. The different parameters of the problem (pressure, velocity, concen- 
tration...) are integrated on some interval of given length (Fig. 3) to provide 

1 The Damkohler number is a dimensionless number, defined as the ratio between the 
acid consumed and the acid transported by convection. 
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the core-scale properties, whose definitions have been given in Part I [2] and 
reminded in the introduction. 

From the resulting porosity, pressure, velocity, and concentration fields, two 
additionnal core-scale parameters are calculated. From the pressure and veloc- 
ity fields, we calculate the head loss between two sections of the core. From 
the concentration field, we calculate the overall dissolution rate for a volume 
included between two cross sections. From these parameters it is possible to de- 
termine the core-scale permeability and the mass exchange coefficient. It must 
be emphasized that their definition is general, and does not depend on a particular 
model. 




: integration length 



L 

Figure 3 - Core-scale approach of the problem. 



Core-scale averaged permeability: 



dp;- 



Core-scale averaged mass transfer coefficient: 



L AP V ' JV, 



These different parameters can be calculated as a function of time for the var- 
ious performed numerical experiments. What is the impact of the dissolution 
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history? In fact, as mentioned in the introduction, a direct relationship between 
the different macroscopic quantities may be historical, due to the coupling be- 
tween the flow velocity and the dissolution process. It must be emphasized 
that this question is of a general interest, and similar problems were found in 
geochemistry, or in dealing with dendritic mushy zones [7, 3]. This problem 
has been solved approximately at the Darcy-scale, where we have already de- 
termined the relations linking the different macroscopic properties, such as the 
permeability versus the porosity or the mass transfer coefficient as a function of 
the cell Peclet number and the porosity. We discuss below the application of the 
same ideas at the upper scale, by focusing particularly on the relations K* - <f> m 
and a* = f (Pe,<j> m , ...). 

3. 1 Macroscopic permeability relationship 

We begin by representing the evolution of the core-scale permeability K* as a 
function of the fluid fraction <p m , or the core-scale porosity ej, for the disso- 
lution regimes: conical wormhole, wormholing and ramified wormhole. The 
reason for keeping the two parameters <j> m and in the analysis have been ex- 
plained in Part I [2]. In the case of local equilibrium, the Darcy-scale porosity 
remains constant in the non-dissolved areas. Thus, the single parameter charac- 
teristic of dissolution is <p m . On the contrary, in the non-equilibrium cases, the 
Darcy-scale porosity varies, and the two parameters (p m and e* p characterize the 
dissolution process, in a somehow independent manner. 

For the two limit cases of dissolution, compact and uniform regime, the prob- 
lem simplifies to a ID case and the dissolution front velocity is constant. It is 
remarkable to see that local equilibrium as well as non-equilibrium dissolution 
can lead to a stationary front displacement. Given the displacement velocity of 
the front is constant, the K* - <f> m relation does not depend on the time evolution 
although it is still difficult to express easily this relationship, especially under 
local non-equilibrium conditions because of the Darcy-scale porosity gradients. 

For the three other cases named above, the question remains open. We tried to 
determine a relation between core-scale permeability and porosity in a similar 
way to the Darcy-scale problem. 
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3.1.1 Compact regime 

For this regime, we have a sharp dissolution front with a displacement velocity 
equal to 

V, = 7r M^ (13) 
— + C„ 

where Vo and Co are respectively the injection velocity and injection concen- 
tration. 

The core-scale permeability field K* is directly linked to the core-scale volume 
fraction 4> m , and their values depend on the position of the cross section versus 
the dissolution front. We can write the relationship K * - (p m , for 0 < (f> m < I, 
under the form 

K* = - (14) 

(1 - (p m ) K fi U i d +(p m K 

where K fi ui j represents the equivalent permeability of the fluid region. 



3.1.2 Conical regime 

Although the effective coefficients A"*(/) and <f) m (t) are a function of time, 
their evolutions, nevertheless, are not independent. Is it possible to obtain a 
relationship between these two parameters which is independent of the dissolu- 
tion history or does it stay a function of time? 
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Figure 4 - Schematic representation of tube dissolution. 



To answer this question, let us focus on a simplified case of our problem. We 
consider a wormhole which propagates into a porous medium of initial perme- 
ability A'o and of width L. The studied part of the domain is far enough from 
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the boundaries so that the flow may be considered steady (see Fig. 4). If we 
represent the wormhole by a tube of width b, the effective permeability K* can 
be estimated analytically. The equivalent permeability for the fluid zone is on 
the order of b 2 /\2 (in 2D) which leads, for the permeability K*,to 



K* = (\-<t> m )K 0 + <t> m b 2 /\2 

* ^y- for 0 OT > 0 

and the channel width b can be also expressed as a function of the fluid fraction, 
i.e. b = <p m L. We obtain: 

K* = ^ (16) 

This relationship gives us an idea about the K* - <f> m correlation existing in 
our problem. 

However, we are generally far from these theoretical conditions. We have 
represented in Fig. 5 the macroscopic permeability as a function of the fluid 
fraction at Q = lcm 3 .h _I for different dissolution times. The integration interval 
is fixed at /,„, = 0.75 cm. It is clear from these results that the permeability 
evolution remains a function of the time dissolution, except for small volume 
fractions <p m , and we have 

AT* = /(tfW) (17) 

This result can be physically explained from the dissolution figures. The 
closest we are from the inlet of the domain (i.e., the bigger the fluid fraction 
is), the less the wormhole can be represented by a tube in the studied section 
and, consequently, Eq. ( 16) is no longer verified. 

3.1.3 Wormholing regime 

On the contrary, for the wormholing regime, the obtained results are closer 
to those predicted by the simplified theoretical model developed as indicated 
above. In fact, in this case, the wormhole geometry is similar to a tube. We have 
represented in Fig. 6, the K* — <j> m relation for different injection rates, different 
values of the mass transfer coefficient a and at different dissolution times. The 
interval length l int is also equal to 0.75 cm. We see that, for this dissolution 
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• Q=1 cm 3 .h 1 - 4=0.5 - C 0 =545 kg.m 3 - 1=16200 s 
Theoretical solution for the capillary tube 
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Fluid fraction <J> m 

Figure 5 - Macroscopic permeability as a function of the fluid fraction (conical regime). 
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Figure 6 - Macroscopic permeability as a function of the fluid fraction (wormholing 
regime). 
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regime, a macroscopic relation quasi independent of the time evolution exists 
and can be expressed as 

K* = f(<f> m ) (18) 

The influence of the integration length has also been investigated. Figure 7 
shows the K* — <f> m relationship, for various values of /,„,. This influence is neg- 
ligible, as long as the interval length remains strictly smaller than the wormhole 
length L w . 
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Figure 7 - Effect of the variation of the integration interval length (wormholing regime). 



3.1.4 Ramified regime 

For the ramified regime, at last, Fig. 8 shows the permeability evolution as a 
function of the fluid fraction for different injection rates, different values of 
the mass transfer coefficient, and different times. It is obvious that a K* — 4> m 
relation does not seem to be a very accurate representation. This is due to 
the impact of the spreading of the dissolution front. In fact, there are now 
three distinct zones in the domain: a fluid zone, a zone corresponding to the 
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initial porous medium, and a transient zone with a porosity gradient due to 
the conditions of local non-equilibrium. The macroscopic permeability K* can 
increase although the fluid fraction value is still 0 (cf. triangles A in Fig. 8). 
Therefore, (p m does not seem to be convenient to correctly represent the geo- 
metry of the domain. While our above discussion suggests that (j> m and e*p 
should be used independently, the results interpreted in terms of e*p only lies 
close to a single curve, Fig. 9. If confirmed, the core-scale relation could be 
written as follows 

K* = f(ef) (19) 



L = 0.75 cm 
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Figure 8 - Macroscopic permeability as a function of the fluid fraction (ramified regime). 

The influence of the length l int on the calculation of K* has also been verified 
and the simulations have confirmed that it was negligible provided it is strictly 
inferior to the wormhole length. 

3.2 Macroscopic mass transfer relationship 

The core-scale mass transfer coefficient a* seems more difficult to estimate. 
Based on the relation at the Darcy-scale a = f (Pe ce u, ep) and from the closure 
problem form, we start by assuming that the coefficient a* depends on the Peclet 

Comp. Appl. Math., Vol. 25, N. 1, 2006 



Cop 



68 



CORE-SCALE DESCRIPTION OF POROUS MEDIA DISSOLUTION 



* 



3 

03 

s 



10 1 



10' 



l„ = 0.75 cm 






• 






• 






AQ=211.5cm\h' 


-A=0.2 


- C 0 =545 kg.m' 3 - 


t=72 s (breakthrough) 


♦ Q=211.5 cm J .h'' 


- A=0.2 


- C 0 =545 kg.m 3 - 


t=63s 


XQ=211.5cm\h' 


-A=0.2 


- C 0 =545 kg.m 3 - 


t=54 s 


■ Q=211.5cm 3 .h' 


-A=0.2 


- C 0 =545 kg.m" 3 - 


t=27s 


• Q=211.5 cm 3 .h ' 


-A=0.5 


- C 0 =326 kg.m 3 - 


t=63 s 


• Q=317.2cm*.rV' 


-A=0.5 


- C 0 =326 kg.m 3 - 


t=30 s 



0.4 0.5 0.6 0.7 0.8 

Macroscopic porosity e.* 



0.9 



Figure 9 - Macroscopic permeability as a function of the macroscopic porosity (ramified 
regime). 

number Pe at the core-scale and the fluid fraction <p m (or the core-scale porosity 
e*p for the local non-equilibrium dissolution regime). Moreover, we assume here 
that the dependence on the Peclet number is small compared with the depend- 
ence versus <p m (or e*p). We will come back later on this assumption. 

We have represented in Fig. 10 the evolution of the core-scale mass transfer 
coefficient a* into the core as a function of the core-scale porosity e*p for dif- 
ferent dissolution times and different injection rates within the conical, worm- 
holing and ramified regimes. This figure shows that the relationship a* — e*p is 
not as evident as the relationship K* — e*p, although the plot is here in semi-log 
scale. Even if a relationship seems to emerge in the case of the conical regime, 
we can remark the same phenomenon of instationarity observed for the K* — ep* 
relationship. At e*p = 0.2, initial porosity of the medium, we recover the mass 
transfer coefficient value at the Darcy-scale. When the fluid fraction (resp. the 
core-scale porosity) increases, the core-scale mass transfer coefficient decreases 
to pass by a plateau before going to 0 when ej — ► 1 . The value depends on 
the dissolution regime. In fact, if a is proportional to \/lp where represents 
the characteristic length at the pore-scale (here = 300 /xm), a* is propor- 
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Figure 10 - Macroscopic mass transfer coefficient as a function of the macroscopic 
porosity e^. 



tional to \jl 2 m where l m is the wormhole width. Therefore, as the wormhole 
width decreases from the conical regime to the wormholing regime, and from 
the wormholing regime to the ramified regime, although it is more difficult to 
define the limit of a wormhole in local non-equilibrium dissolution, the value of 
the core-scale parameter a* decreases in the same way. Thus, when the macro- 
scopic porosity s*p increases, the core-scale mass transfer coefficient decreases 
of a factor (lp/l m ) versus its initial value at el = 0.2. 

The uncertainty about the relationship for a* raises a major problem: like at 
the lower-scale, several numerical simulations confirmed that the propagation 
velocity of the wormholes predicted by our core-scale model is strongly sen- 
sitive to the value of this coefficient, especially in the zone at the end of the 
wormhole (e^ = 0.3). Based on the assumption that this uncertainty is linked 
with the multitude of channels created at the short times at the core inlet, we 
decided to eliminate from the porosity and concentration fields the wormholes 
that do no longer develop, and to only consider the influence of the dominant 
wormhole (this will be referred to as the dominant wormhole assumption). The 
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Figure 1 1 - (a) porosity field and (b) a* — correlation in wormholing regime before 
the dominant wormhole assumption; (c) porosity field and (d) or* — s*^ correlation in 
wormholing regime after the dominant wormhole assumption. 



interest of using this assumption is confirmed by Fig. 1 1 . The variation of the 
integration interval length has always a negligible influence on the correlation. 
An interesting parallel can be made between this assumption and the quasi- 
stationarity condition often used for solving the closure problems. With this 
last assumption, non-local effects in time and space are eliminated. In the same 
way, neglecting these wormholes allows to neglect the spatio-temporal varia- 
tions of the deviations, and, consequently, the theory is valid only for some long 
times and far enough from the core inlet. 

A zoom around the zone where the mass transfer coefficient strongly de- 
creases, and which contains the most important part of the dissolution physics, 
is represented in Fig. 12. This shows a non negligible dependence of the a* 
coefficient versus the injection velocity, i.e., the Peclet number. It is therefore 
necessary to take into account both the fluid fraction (or the porosity and 
the injection Peclet number for the development of the correlations for the core- 
scale mass transfer coefficient. The log-normal representation in Fig. 12 sug- 
gests the use of correlations for a* with an exponential dependence with e*^ 
the magnitude varying as a function of the injection Peclet number. 



Comp. Appl. Math., Vol. 25, N. 1, 2006 



FABRICE GOLFIER et al. 



71 




0.07 0.08 0.09 



0.1 



♦ Q=21.1 cm 3 .h 1 - C 0 =545 kg.m 3 - 1=450 s (l int =0.75) 
□ Q=50.3 cm 3 .h" 1 - C 0 =545 kg.m 3 - 1=152 s (U=0.375) 

■ Q=50.3 cm 3 .h' 1 - C 0 =545 kg.m* 3 - 1=152 s (1^=0.75) 
o Q=50.3 cm 3 .h" 1 - C 0 =545 kg.m' 3 - 1=76 s (l int =0.375) 

■ Q=50.3 cmV - C 0 =545 kg.m 3 - 1=76 s (l int =0.75) 
o Q=211.5 cmV - C 0 =326 kg.m' 3 - 1=63 s (l int =0.375> 

• Q=211.5 cm 3 .h 1 - C n =326 kgjn" 3 - 1»63 8 (U s 0 75) 



Corr6lat|QP? 

Q=21.1cm 3 .h' 

Q=50.3 cm 3 .h 1 

... Q=211.5cm 3 .h 



Fluid fraction (J>. 



Figure 12 - Expression of correlations or* = / (Pe, <p m ) 



As a conclusion, several important results can be extracted from the develop- 
ment of these correlations: 

• In the conical regime, effective coefficients cannot be uncoupled from 
the dissolution history. 

• Simple correlations for the mass transfer coefficient remains difficult to 
be evaluated numerically. This represents an important problem for the 
development of a macroscopic model at the core-scale, while it is not 
clear at this point what is the impact of errors made in determining these 
correlations. 

• The obtained correlations do not depend on the length of the integration 
interval /,„,, within the range of parameters explored. 

4 Core-scale ID model 

From the system of equations obtained at the core-scale and the mass trans- 
fer and permeability correlations provided by the numerical simulations at the 
Darcy-scale, two ID core-scale models can be used to verify the validity of 
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these representations. The one-medium or one-equation local non-equilibrium 
model corresponds to the coupling between the flows equations, Eqs. (3) and 
(4), and the transport and dissolution equations, Eqs. (8) and (9). The devel- 
opment of the two-medium or two-equation local equilibrium model, based on 
the reaction-transport equations presented previously, Eqs. (6) and (7), needs 
the introduction of an additional assumption. In fact, this model requires the 
knowledge of the regional average velocities and V*, i.e., solving of the 
regional core-scale Darcy equations written as 



v; = o (20) 



V; = - -K; • (V^ - Pp g) in the m - phase (21) 



and 



V • VJ = 0 (22) 

VJ = - -KJ • (VP; - p p g) in the n - phase (23) 
M 

where and K* represent the regional effective permeability tensors [8]. 
The two average velocities are linked with the following relation: 

= 0arU; + (l-0 w )U; 

The intrinsic average velocity in the porous medium U* being negligible ver- 
sus the velocity in the fluid region U^, we can assume 

v; * v; (25) 

which allows us to use Eqs. (3) and (4) for the flow part. This approximation 
raises some difficulties only for some small values of <\> w , when the wormhole 
is not yet developed. 

In this section, we test the capacity of these models to describe the core-scale 
dissolution physics. For this purpose, we compare the predictions given by the 
1 D averaged models with the simulations performed at the lower scale. Since 
the large-scale dispersion tensor DJ*, which appears in the transport equation, 
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201 x 101 nodes 


A =0.5 


p a =2160 kg.m~ 3 


/x = 1.10~ 3 Pa.s 


e*=0.2 (or0* = 10" 4 ) 


P = 1.62 


D = 2.10~ 9 m.s~ ] 


K* = 1.5- 10- 11 



Table 1 - Numerical data for constant flow rate injection. 



has not been calculated, it has been taken equal to <p m D (or e^D for local non- 
equilibrium dissolution) in our simulations. This assumption has a small effect 
on the propagation time because we have focused here on the wormholing and 
ramified regime, for which the Peclet number is relatively important. The Darcy- 
scale simulations which are used as reference, are obtained from the model 
described in Golfier et al. [1] and described briefly in Section 2. The local 
characteristics of the porous medium and its core-scale properties, including the 
initial values of e J (or <j) m ) and K * are given in Table 1 . The correlations used for 
the evolution of the mass transfer coefficient a* and the effective permeability 
K* are obtained as previously described. 
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Figure 1 3 - Comparison of the Darcy-scale model and the core-scale models for Q = 2 1 
cm 3 .h~' and Co = 545 kg.m -3 . 



The comparison of our results for both models with Darcy-scale simulations 
are illustrated in Fig. 13 and 14 for Q = 21 cm 3 .h _l - Co = 545 kg.m -3 and 
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Figure 14 - Comparison of the Darcy-scale model and the core-scale models for Q = 50 
cm 3 .h _1 and Co = 326 kg.m -3 . 

Q = 50 cm 3 .h -1 - Co = 326 kg.m -3 , respectively. It must be emphasized 
here that the choice of the porosity threshold used to determine the tip of the 
wormhole is a decisive criterion: it has been fixed from the observation of the 
porosity fields at e* p = 0.25 (<f> m = 10~ 2 for the two-medium model). In the 
case of the two-equation model, the initial value of (f) m at t = 0 s, theoretically 
equal to 0, is fixed at 10 -4 for numerical reasons. The one-equation model does 
not succeed to correctly reproduce the propagation velocity of the wormhole 
and overestimates the breakthrough time in both regimes. The wormhole prop- 
agation, indeed, is represented in the one-equation model by an increase of Ep 
which is a property averaged on all the domain under consideration. For in- 
stance, a uniform dissolution of the porous medium and a conical wormhole 
well-developed could lead to the same value of e* p in certain conditions. As 
a consequence, the macroscopic porosity and hence the one-equation model, 
are not suitable to describe the local behaviour of the wormhole developement 
within a porous medium. On the contrary, results are more promising for the 
two-equation model. The propagation velocity is correctly predicted (same slope 
coefficient), except for the short times where the propagation time is underes- 
timated. This difference can be explained by the fact that the phenomenon of 
wormhole competition is not taken into account in the calculation of a* (dom- 
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inant wormhole assumption). This emphasizes the difficulty to express this 
coefficient independently of the history. The comparison of the fluid fraction 
fields obtained, either directly from the core-scale model, or by spatial integra- 
tion of the Darcy-scale porosity fields, confirms this assumption. One can see in 
Fig. 1 5, that the (p m value is correctly predicted, except for the inlet of the domain 
where it is strongly underestimated because the core-scale model does not take 
into account the multitude of channels which propagate over a few centimeters. 




Figure 15 - Fluid fraction values predicted by the Darcy-scale model and the core-scale 
2-equation model for Q = 50 cm 3 .h -1 and Co = 326 kg.m~ 3 at / = 152 s. 

These first results suggest that it is impossible to completely eliminate the 
non-local effects and that a unique correlation for the mass transfer coefficient, 
independent of x and t, does not allow to perfectly reproduce the dissolution 
phenomenon at the core-scale. In the framework of the proposed core-scale 
models, several routes can be considered to model these mechanisms and correct 
the inaccuracy induced by such a representation: (i) taking into account the time 
effects by the introduction of a different a* -correlation at short times such as 
a * = ^compact — constant for t < t 0 , (ii) taking into account the spatial effects 
by the introduction of a different a* -correlation at the inlet of the domain such as 
a * = a compact = constant for x ^ xq. We have chosen here, from the observa- 
tion of the fraction fluid values, a "mixed approach" by using a different corre- 
lation for a* where the porous medium at the inlet is not dissolved enough, i.e., 
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a * = Compact = constant for <j>™ let < (f> 0 . Then, the a* ompacl variable becomes 
an additional parameter of the model that may be fitted to "experimental" data. 
The results of the introduction of this inlet compact dissolution model are illus- 
trated by the curve with the symbols (•) in Fig. 13 and 14 and the impact on the 
fluid fraction field is presented in Fig. 16. 




Figure 16 - Fluid fraction values predicted by the core-scale 2-equation model with 
compact front for Q = 50 cm 3 .h -1 and Co = 326 kg.m -3 at / = 152 s. 

These different results confirm that the description of the studied system by a 
two-medium domain may be an appropriate way for describing the dissolution 
phenomenon at the core-scale. 



5 Conclusion 

We have determined in this paper the correlations used for the effective coef- 
ficients which appear in the two core-scale transport-reaction models developed 
in Part I [2]. In both cases, the correlations have been directly extracted from 
the study of Darcy-scale numerical simulations and not from the solution of 
"closure problems". 

Several questions are left open concerning both the form of the macroscopic 
equations and the correlations used for the effective coefficients. Nevertheless, 
the results suggest that a two-medium model is a better candidate to describe 
the physic behaviour of the dissolving system. Furthermore, it seems that the 
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effective coefficients cannot completely be uncoupled from the non-local effects 
and the expressions of their correlations remains non-local functions, i.e., func- 
tions of x and t. A good estimation of the mass transfer coefficient a* allows 
however to correctly reproduce the wormhole propagation into the domain. 

What are the implications of these results for larger "averaging" volume, i.e., 
field scale problems? The introduction of a supplementary upscaling to de- 
scribe the wormholing phenomenon at the large-scale seems difficult. The use 
of our core-scale model as a near well-bore simulator coupled with a classical 
field model seems preferable. In fact, the results of our model could be used to 
calculate the skin effect factor in order to take into account the real conditions 
around the well within the field model. Therefore, we investigate at this time 
the ability of core-scale models to correctly reproduce the experimental results 
being based on accessible experimental data only, i.e., porosity and fluid fraction 
fields which can be obtained from tomography methods. 
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Abstract. The versatility of a robot to perform a task is limited principally by the flexibility 
of its end-effector. In the last years, research has been focused on the development of a hand with 
several fingers since these devices are capable of manipulating and grasping objects of different 
forms. A dexterous manipulation system, composed of a robot hand with several fingers and 
an object that will be held or manipulated, could be modeled as a set of rigid bodies in contact. 
The dynamics of several rigid bodies in contact tries to predict the accelerations and forces at the 
contact points of the set of rigid bodies with Coulomb friction. The calculation of such forces 
allows us to determine if the contact is maintained or disappears and to plan a determined action. 
The equations that describe the problem form a system of differential algebraic equations. In this 
contribution the problem is reformulated as a mixed nonlinear complementarity problem (MNCP). 
Then, an optimization problem with box constraints associated to the MNCP is presented using an 
adequate merit function. Conditions about the equivalence between the problems are established. 
Finally, the optimization problem is solved using a robust and efficient algorithm. Encouraging 
numerical results are reported. 

Mathematical subject classification: 49M37, 65C20, 90C30, 90C33. 

Key words: linear complementarity problem, box constrained minimization, multi-rigid-body 
contact problem. 



1 Introduction 

Complementarity problems are of great importance in engineering applications 
because they are associated to the notion of equilibrium systems [4]. In the last 
decades, different types of complementarity problems have been studied: linear, 
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nonlinear, mixed, and efficient and robust algorithms have been developed to 
solve each of them. 

This kind of problems arise in many engineering applications [4], such as 
mechanical contact problems, structural mechanics problems, structural design 
problems, nonlinear obstacle problems, elastohydrodynamics lubrication prob- 
lems, traffic equilibrium problems and optimum control problems, to mention 
some. 

This work is organized as follows: section 2 presents generalities of the mixed 
nonlinear complementarity problem, and section 3 presents both the method 
for solving the MNCP via optimization, and existence and uniqueness results. 
Section 4, considers a dynamic model for several rigid bodies in contact, the 
equations that represent it [8, 9] and from these, the formulation of MNCP. 
In section 5 the numerical results for three-fingered hand manipulation system 
holding an object are reported. Finally, conclusions and future possibilities to 
continue on this research line are presented. 

2 Mixed nonlinear complementarity problem 

Let J : R" — ► R n be a vectorial function and £ and C be index sets that define a 
partition of { 1 , 2 n). The problem of finding a vector jc e £2 c R n such that 

f L (x) = 0, x c >0, fc(x)>0, x' c f c (x) = 0, (1) 

is called the mixed nonlinear complementarity problem. The variables x c are 
the complementarity variables and X£ are the free variables that do not satisfy 
the complementarity and non negativity conditions. 
The set Q is defined by 

ft = {x e R" : a < x < b}, 

where a and b are n -dimensional vectors with a t € [-co, oo) and € (a, , oo]. 
If Xi is such that / e L then a, , = -oo and b t = oo, while if / € C then a { = 0 
y hi = oo. 

If the set L = 0 and Q = R n + = {x e R" : x t > 0, i = 1 n) then 

we recover the nonlinear complementarity problem (NCP), if the function f(x) 
is also affine then it is a linear complementarity problem (LCP). If the set C = 0 
it results in a nonlinear system of equations. 

A change in notation is considered, introducing a partition of x, the vectors 
u € R p , v £ R m and a partition of J the vectorial functions F : R" R p , 
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and G : R" — ► R m . Slack variables z e R p are introduced, and the mixed 
complementarity problem is expressed as follows, 

«, z > 0, z- F(u, v) = 0, u'z = 0, G(u, v) = 0, (2) 

where u is the vector of complementarity variables and v is the vector of free 
variables, so that the dimension of the complementarity problem is m + 2p. 

The mixed nonlinear complementarity problem can be considered as a non- 
linear complementarity problem with an extra simultaneous nonlinear equa- 
tions system. 

3 Solving of the complementarity problem via optimization 

In order to reformulate the mixed complementarity problem as an optimization 
problem, it is necessary to introduce a merit function / : R" x R' M x R p -> R. 
The optimization problem that is solved is 



The advantage of reformulating the MNCP using the aforementioned strategy 
is the possibility of applying efficient algorithms and theoretical results known 
for optimization problems. In general, the algorithms used to solve optimization 
problems obtain stationary points, that is, points that only satisfy the first order 
optimality conditions. The convergence to a global minimum can be guaranteed 
if the merit function is convex. 

Optimization problems have been extensively studied in recent years and many 
efficient algorithms exist to find the solutions. To obtain a solution to the mixed 
nonlinear complementarity problem (3) with the merit function, 



Andreani, Friedlander, Mello and Santos [1] establish conditions that allow 
us to assure that stationary points of (3) with the merit function (4) are also 
global minimizers and solutions of the MNCP. 

4 Mathematical model of the contact problem 

To work with contact problems with friction, a model of surfaces in contact and 
bodies in movement is considered. 
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u > 0, z > 0. 



(3) 



/(«. v. Z ) = II F{u, v) - z || 2 + || G(«, v) || 2 + (u'z f , 



(4) 



Copyrighted material 



82 



A MIXED NONLINEAR COMPLEMENTARITY TECHNIQUE 



The dynamic system is formed by a determined number of passive bodies 
called objects and a given number of active bodies called manipulators. The 
first ones move in response to external and contact forces. The hypotheses of 
the proposed model [7, 8] are: 

1. The bodies are rigid. There are no constraints on the shape of the bodies. 

2. The normal direction at each contact point is well-defined, this means that 
the surface of the passive body has a tangent plane on each contact point. 

3. There exists friction at each contact point. A friction model and a friction 
law must be selected. 

4. The manipulator is formed by links and joints; the manipulator joints do 
not form closed loops. Each joint has only one degree of freedom. All 
allowed contacts are unilateral. 

5. All links are connected to a grounded link (reference coordinate). 

The equations and constraints that govern the model are: the Newton-Euler 
movement equations, for the objects and manipulators, the unilateral and bi- 
lateral kinematic constraints due to the contact, the dynamic conditions of the 
contact, that is the non interpenetration principle and non tensile forces, and 
the friction law of the contact. 

To develop the equations that govern the dynamics of the set of rigid bodies 
in contact, it is considered that in a certain instant t 0 certain number of objects, 
n obj , are in contact among them and with a certain number of manipulators, 
n m an, being n c the total number of contacts. 

Let Cij denote the contact force acting on the body i through contact j expressed 
in the contact frame C 7 , with normal, tangential and orthogonal components 
icij) n , (Cij), and (c, y ) 0 . The sum of the forces acting on the body i is equal to 
the mass times the acceleration of its mass center and it is expressed as 



where % is the set of indices of the points of contact of the body i t g obj j e R 6 
is the external generalized force (expressed in the fixed reference coordinates of 
the body B, ), which acts upon it, h obj<i e R 6 is the term of movement quantity, 
q ob jj € R 6 is the generalized linear and angular acceleration of the mass center 
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of the object, M obj j g R 6x6 is the symmetric positive definite mass matrix of 
the object i partitioned as: 

= [ M o' l } ■ 

with Mi j = m obj j / 3 , / 3 being the identity matrix, and Af 2 ,, the bodie's inertial 
tensor / expressed in the base B,. The matrix W u e R 6x3 transforms the contact 
forces Cij into a generalized force, F u . The wrench matrices Wjj contain all the 
geometric information about contact j. 

The previous equations extended to all objects, express the contact forces by 
components, and the following matrix equation is obtained: 

W n c n + W, c, + W 0 c () + g obj + h obj = M obj q obj , (7) 

with c a <e R"<\ W a e R 6 ^ x,lf ,ff e {n,t,o}, the vectors g obj , h obj , q obj e 
R 6 "^, and M obj e R 6 6 is 

M obj = diag{M objA , M obj . 2 M objMobj }. (8) 

The dynamics equations for the manipulators are obtained likewise: 

r - (7,; c n + // C t + J l 0 c fl + g man + h man ) = M man 6 man , (9) 

where, if we call n 0 the number of manipulator joints, G man e R n " is the velo- 
city of the joints, r e R n ° is the vector of joint effort, J n , J, and J a e R w « x "' 
are the Jacobian matrices which determine the effect of the a-th component 
of the force (c,y)« over the component a of the vector of joint effort r, with 
a € {n, r, o}, M man e R"" x "" is the symmetric positive definite inertia matrix, 
gman Wman) e R"° is the vector of torques induced by gravity, h man (B man , 0 man ) e 
R"° represents the vector of torques and forces induced by Coriolis and centri- 
petal accelerations and 9 man e R n " is the vector of joint accelerations. 
The motion of each object in contact is subject to kinematic constraints: 

Va = W^ohj - J a 0man, Of <E {//,/,(?}, (10) 

<*a = W^ijohj - J a d man + W^q obj - j a Q man , OL G {«, t, o}. (11) 

where v a €. R" r is the linear velocity and a a e R" c is the linear acceleration in 
the direction a, expressed in the contact frame C,-, q ob] and B man are the vectors 
of the velocities of the bodies and manipulators. 
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In each contact point j the normal component of the velocity is zero, 

v jn =0, ; = l,...,n c . (12) 

Among the n c contact points, « R denotes the number of rolling contacts and n$ 
denotes the number of sliding contacts, so that n c = ns + n-R. S and H denote 
the sets of respective indices. 

A rolling contact point is characterized by the fact that both the tangential and 
orthogonal components of the velocity at such point are zero, that is 

v Jt = Vjo = Q, /eft, (13) 

and a sliding contact point is characterized by the fact that the tangential or 
orthogonal components of the velocity at such point is different from zero, 

vjt^Q 6 tv 7*0, jeS. (14) 

To avoid interpenetration, the motion of the object is subject to the constraint 
known as non penetration principle ofSignorini, 

a n > 0. (15) 

The forces on the contact point must have non negative normal components, 
this means that they are not tensile, 

Cj n > 0, j = 1, ...,n f . (16) 

For any contact point j, if (a„)j = 0 the contact is maintained and (c„) y > 0, 
whereas if (a„)j > 0 the contact disappears (breaking) and (c n ), = 0. 

Considering (15) and (16) relating the accelerations and the forces of the con- 
tact points, the complementarity constraint is established, 

a' n c n =0. (17) 

To grasp an object, friction forces are required and so a friction model must 
be established. The choice of the friction model determines the form of the 
matrices W a which appear in the equation (7) and the matrices J a in equation 
(9) with a € {n, t, o}. The law used here is the Coulomb friction law which 
establishes that the contact force on the point j falls within or on the boundary 
of the corresponding friction cone, represented as follows: 

Cj t + c 2 jo < n 2 jC 2 jn , j = 1, n c , (18) 
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where fij is the friction coefficient on the contact point j. 

If the contact is sliding, the contact force must lie on the border of the friction 
cone with its friction component directly opposed to the sliding velocity and 
satisfy 

(19) 



Cj n Vj a + CjaJtf, + vj g = 0, j 6 S, Of € {*, o). 



If the contact is rolling, the contact force can have some direction and magni- 
tude, within the cone defined by (16) and (18) and satisfy 



lijCj n a ja + Cj ay ja] t + a 2 jo = 0, j € !R, a € {/, o). 



(20) 



Equation (19) differs significatively from equation (20), since the former is 
linear in the unknowns cj„ and Cj a because the velocities are data, whereas the 
latter is nonlinear because the accelerations are unknowns. 

Considering that in the initial instant to q ( ,i,j, q 0 hj, gobj, h on j, M on j are 
known for the objects; 0 man , G man , r, g ma „, h man , M man for the manipulators; 
the torque matrices W and the Jacobian matrices J man , q 0 bj and 6 man satisfying 
the kinematic constraints (10), (12) and (13), to solve the dynamics of several 
bodies in contact means to determine q 0 ^ i, 9 man ,c n , c t , c Q , a n , a { , a Q , which 
satisfy the equations (7), (9), (1 1), (15), (17), (19) and (20). 

Since the algebraic-differential equation system that results is complex to 
solve, it is convenient to reformulate it as a mixed nonlinear complementar- 
ity problem. Thus, a discretization of the independent variable is considered, 
assuming that in the instant to an initial configuration, i.e. a contact and the 
position and velocities of all the bodies is known. Therefore, to establish the 
configuration in an instant to + A/, a complementarity problem is solved. This 
implies analyzing whether the contact is maintained or not, and calculating the 
positions and velocities of the new configuration. 

Taking into account that the inverses of M„hj and M man exist, q„bj and 9 man can 
be obtained from the equations (7) and (9). Replacing q () bj and 0 man in equation 
(11) and defining a a e W c , with a e [n, t, o), a linear equation system is 
obtained 



(21) 
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" c n 




' b„ ~ 


a, 


= A 
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A f> n 


A 0 t 


A ( ,o 



= J'MJ, 



(22) 
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where the A ay are submatrices of A which consist of rows and columns corre- 
sponding to the directions a and y respectively, with a, y e {n,t,o}. 
From M obj , M man , W n ,W t , W Q , //, y J l 0 , the matrices 



M = 



(M obj ) 
0 



-l 



(M, 



o ] _ r w n w t w 0 ] 

J ' 3 - V 4' * A 1 J ' (23> 



and the vector 



= j'[ > l+J'^ **>t*«* 1. 

L "man J |_ 6/w<jw > "wia/i T J 



(24) 



are denned. 

Equation (19) is solved to obtain c,, and c 7< , with j e S and replace them in 
(21). The acceleration vectors and contact forces are partitioned for the contacts 
R and S. Since there are no additional constraints for the tangential and orthog- 
onal components of the accelerations on the sliding contact points, the equations 
of the system (21) which define these components can be set apart without af- 
fecting the solution of the complementarity problem 



(25) 
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(26) 



~(A n n)sS~ 




~(Ann)ss~ 




~(A n t)sS~ 




~(A no )ss~ 














(Ano)'RS 












Vts - 


{A to )RS 






_(A 0 n)ns_ 




_( A ot)ns. 




_{Aoo)'RS_ 



Vos, (27) 



V aS = diag 




My v ja 



, j € S, 



+ v 



jo, 



a € {t, o}. 



(28) 
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Equations (15-17), (18-20), (25-28), constitute a mixed nonlinear comple- 
mentarity problem as defined in (1) associated to the 3D dynamics of several 
rigid bodies in contact with Coulomb friction. 

To express the mixed nonlinear complementarity problem with the notation 
used in (2), the slack variables sj are introduced by 

and the vector of complementarity variables X whose components are given by 

(30) 



Vector u e R ns+2n *- corresponds to the complementarity variables, vector 
v e R 4 " R corresponds to the free variables and the vector of slack variables 
z e R n * +2n * are defined 



u = 



CSn 
C'Rn 

X 



> 0, v = 



C'Rt 
C'Ro 

dRo 



Z = 



ClRn 

s 



>o. 



(31) 



Partitioning adequately the system (25) and from (29), the following vectorial 
function F : R"s +6 "k R n s+ 2 "* is obtained, 



/ 



F(u,«) = 



+ 



CSn 
CR„ 
CRt 
CRo 

S(CR„, C'R t , o Ro ) 



\ 



bRn 



(32) 



with 



L (A nn y RS (A wn ) KR J ' 2 " [ 



(A n ,)RR (A no 



no)sR j 
no)RR J 



The vectorial function G : JR"** 6 "* R 4 "« is expressed as follows 



G(u,v) = 



G\(u t v) 
G 2 (u,v) 
_ G 2 (u,v) _ 



(33) 



where 



(G,(w, v))j = HjC jn a jt + Cjtkj, j € R, 



(34) 
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(G 2 (m, v))j = fijCj n a jo + Cj 0 \j, j e R, 



(35) 



* - [ a 
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C'Rn 
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r b nt ] 



)RR (A to )RR 
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] 



(36) 



5 Numerical results 

In this section, we show how the mixed non linear complementarity technique 
presented in the previous sections is used. Let us consider the model of an object 
grasped by a robot hand with three fingers. 

The initial configuration for the proposed example is established fixing the 
location of the coordinate systems with origins in the mass center of the body 
(B), the base of the j-th robot (S ; ), the j-th point of contact (C y ) and the palm 
of the hand (P). The origin of the reference coordinate system is located in 
the palm of the robot hand. Figure 1 presents a bidimensional scheme for the 
given configuration. 




Figure 1 - Three fingers holding a sphere. 
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Using algorithm 1, the system configuration of various bodies in contact in the 
instant t 0 + At from an initial configuration is obtained. 

Algorithm 1. Given p, n< R , n s , and the configuration parameters, 

1 . Define the initial configuration. 

1.1. Calculate the coordinate transformation matrices, 

1 .2. Calculate the vector of positions of the manipulator joints, 9 man 
(Problem of inverse kinetics). 

2. Calculate the matrices for the object W, M„bj. 

3. Calculate the manipulator matrices, J man , M man . 

4. Calculate the matrices of the complementarity problem, M, J, and A. 

5. Solve the complementarity problem via optimization. 

6. Define the new configuration. 

In step 5, the optimization problem can be solved by any appropriate algorithm. 

To establish the initial configuration, the matrices p t and p u € E 3x " ( are 
considered. Each column of p t , /?,(., j) contains the j-th point of contact with 
respect to the system P; while the j-th column of p„ contains the coordinates 
of the 7-th point of contact with respect to the coordinate system S 7 . For the 
proposed configuration, the matrices are, 



Pt = 



r cos{o) r cos(o) —r 
—rsin(o) r sin{o) 0 
a cm a cm a cm 



Ptt = 



-(b - r)cos(cr) -(b-r)cos(o) b — r 
(b-r)sin(a) -(b - r) sin(cr) 0 
acm a cm a, 



cm 



where r is the sphere radius, a cm establishes the location of the sphere mass center 
with respect to the coordinate axes system of the palm P, p cm = (0, 0, a cm ), b 
is the distance between the origin of the coordinates B and the projection on the 
plane z = a cm of the origin of the y-th coordinate system S y , j = 1, 2, 3 and 
o is the angle formed by the positive axis x of the Cartesian coordinate system 
fixed in the mass center of the body B and the segment which joins the origin of 
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B with the points of contact Ci and/or C 2 . As observed in figure 1, the location 
of the contacts Cj and C 2 is symmetric respect to B. 
The rotation matrices for each contact have the following expressions: 

Kp.c, = Qz(-cr) Q y (^y /?p,c 2 = Q*(j) &(-*>. 

where Q a (<p) are the matrices representing elemental rotation around the direc- 
tion a, with an angle (p. 
The friction model is characterized by the matrix, 



1 0 0 

0 1 0 

0 0 1 

0 0 0 

0 0 0 

0 0 0 



(37) 



The matrix W € R 6x9 is expressed as follows, 

W = [ W x B f W 2 B f W 3 B f ], 
with Wj, with y = 1,2,3, defined as, 

Wi = 



(38) 



r Kp.c, o i 



(39) 



The following notation is used: given any vector v e R 3 , v e M 3x3 we define 



v = 



0 

v 3 

-v 2 Vi 



-v 3 
0 



v 2 

-V] 

0 



(40) 



The inertia matrix of the sphere respect to the mass center is M 2 = \mr 2 l 3 , 
where l 3 is the identity matrix of dimension 3. When the fixed reference system 
is not in the mass center of the body a modification of the inertia matrix must be 
performed using the result of the parallel axes theorem [3]. 
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The mass matrix of the object M obj gK 6x6 is, 



Mobj 



I" m/ 3 0 1 
= L 0 <>J 



(41) 



The experiments were performed using a well known robot PUMA (Pro- 
grammable Universal Manipulator for Assembly). For each manipulator the 
Jacobian matrices Jf. and the mass matrices M f, are calculated using the robotic 
toolbox in MATLAB [2]. The Jacobian matrix for all the three fingered hands 
J man e R 9xn ° results 

J\ 0 0 
0 h 0 
0 0 h 



Jman — 



(42) 



with 



where B/ is the matrix of the friction model. 
The inertia matrix of the manipulator M Sl e R"'i *"*> is 



M fj (9) = Y,J t L J (0)ML J JL j (O), 



(43) 



where // and Mi are the Jacobian matrices and the generalized inertia matrices 
of the link Lj. The mass matrix of the three fingers M man G W" xn ° is 



M fl 0 0 
0 M h 0 
0 0 M h J 



(44) 



If it is considered that all the contacts are sliding, the MNCP is reduced to a 
LCP. Solving the LCP means solving the optimization problem defined in (3) 
with the merit function 



2 ■> 

f(c n ,a n ) = || a n - A nn c n - b n || 2 + (a' n c n ) , 
where the vector b n e R 3 is obtained from equation (24), 



(45) 



A nn — A nn A nt V ( A no V Q . 



(46) 
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The matrices V t and V a are defined in (28) and each submatrix A ay has the 
expression 

A ay = Wa(M 0 bj)~ l Wy + J a (M man )- l J Y '. (47) 

A variable change is performed, t e R 6 , t = [c n , a„Y so that the merit function 
is expressed as follows, 



3 / 3 

/(/) = £h>3-E Afyfy - (Mi ] + 



(48) 



The matrix M in the complementarity problem depends on the friction coeffi- 
cients on the contact points. For the proposed example, and using a PUMA560 
manipulator the results obtained are: 





Boxlt 


FunEval 


Quacanlt 


Matxvec 


/ 


0.1 


4 


5 


15 


22 


2.0526 x 10- 23 


0.2 


6 


7 
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Table 1 - Results of the optimization problem for the LCP. 
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Table 2 - Results of the LCP for the sliding contacts. 

As can be observed in Table 2, the contact is mantained in all cases. If it is 
considered that the three contacts are rolling, the problem is formulated as a 
MNCP. The MNCP is solved considering the equivalent optimization problem 
(3), using the merit function proposed in (4). The following variable change is 
performed t = [c„, c,, c 0 , a n , a t , a 0 , s , A.]', with the merit function 

f(t) = sum\ + sum 2 + sum] + 5wm 4 , 
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Table 3 - Results of the optimization problem for the MNCP. 
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Table 4 - Results of the MNCP for the rolling contacts. 



where 

3 2 

sum x = J2[hj + 18) - (/xO) t{j)) 2 + (t(j + 3)) 2 + (t(J + 6)) 2 ]", 

3 2 

sum 2 = J^[nU) t(J) hj + 12) + t(j + 3)7 U + 21)] 

+ [ti(j) hj) hj + 15) + hj + 6)hj + 2i)] 2 , 

3 3 

sum 2 = hj)hj + 9) + J2 t(j + mhj + 21), 
sum 4 = \\z\- Fi(«, v) \\\ + || G 3 (w, v) \\\ . 
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The optimization problem solved is of dimension 24. It has been solved by 
using a robust and efficient algorithm developed by Friedlander and Martinez 
[5, 6]. The results obtained are summarized in tables 3 and 4. 

6 Conclusions and future research 

The dynamics of several rigid bodies in contact with Coulomb friction is pre- 
sented as a MNCP. The complementarity problem is reformulated as a box 
constraints optimization problem. 

It is worth noting that the optimization problem for the mixed case is general. 
Contrary to many methods which are specific to linear or nonlinear complemen- 
tarity problems, this technique can be applied to any of them by only defining 
appropriately the merit function and the box. 

An example of a robot hand with three fingers grasping an object is considered. 
Matrices of smaller dimensions, which respond to a real model, were preferred 
over those with greater dimensions. Functions in MATLAB were developed to 
obtain the system matrices and the optimization problem with box constraints 
was solved using the easy computing program. The robustness of the optimiza- 
tion method we used, allowed us to obtain a solution in most cases. In the future, 
we wish to generate contact problems with a greater number of bodies and points 
of contact, both rolling and sliding, trying to consider the structure of the matrices 
which are generated without losing sight of the physical meaning. Furthermore, 
it is considered of interest to complete the kinematics and dynamics calculation 
of the robot problem for a sequence of times t\, ti, t^... and to verify if the 
obtained results are satisfactory. We plan to continue deepening the study of the 
MNCP in general following the proposed research line. 
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1 Introduction 



In this paper we consider the following convex quadratic program in its classical 
standard (QP) form called the primal 



minimize 

X 



^c T x + ^jc r <2*j , s.t. Ax = b, x ^ 0, 



(1) 



and its dual form 



maximize 

y,x, z 



b T y - ^x T Qx 



, s.t. A T y + z - Qx = c, O 0, (2) 
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PATH-FOLLOWING METHOD FOR QUADRATIC PROGRAMS 



where Q is a given (n x n) matrix, c e W, fteR" 1 and A is a given (m x n) 
matrix. 

Recently, some other important formulations have been used to study con- 
vex quadratic program such as the (second-order) cone optimization problem 
(SOCP). They are currently the most often used for quadratic convex programs. 
See, e.g., [5]. 

Quadratic programming is a special case of nonlinear programming and has a 
wide range of applications. A classical method to solve convex (QP) is the Wolfe 
algorithm which can be considered as a direct extension of the well-known sim- 
plex method, in the worst case its complexity is exponential. Therefore, there is a 
growing interest for developing efficient, robust and polynomial time algorithms 
to solve the problem. Primal-dual path-following methods are among the most 
efficient algorithms of interior point methods (IPM) to solve linear, quadratic 
programming, complementarity problems (CP) and general conic optimization 
problems (COP). These algorithms are, on one hand, of Newton type and thus 
lead to an efficient practical implementation and on another hand they have a 
polynomial time complexity. Recently, many interior point methods have been 
proposed for solving (QP) problems in its standard format or via (SCOP). These 
methods are based on new search directions with best theoretical properties such 
as polynomial time complexity and efficient practical implementations. See, 
e.g., [2, 4, 5, 6, 7]. In this paper, we extend a recent approach developed by 
Darvay [6] for linear programs to (QP) case. In this approach, based on the 
introduction of a monotonic function <p of one real variable. The equations 
x,Zi = it in the nonlinear system which defines the central path (see Wright [8]) 
are replaced by 



Of course, if <p is the identical function we recover the classical path-following 
method. For (p(t) = V? we define a new class of search directions and hence a 
new primal-dual path-following method for approximating the central path. We 
prove that the related short-update algorithm solves (QP) locally and quadrati- 
cally and in a polynomial time. We mention also some related and interesting 
IPM approaches studied recently by Bai et al. [4]. They have also defined a new 
class of search direction for linear optimization by using the so-called kernel 
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function, that is any univariate strictly convex function k(t) that is minimal at 
/ = 1 and such that k(\) = 0. However, by defining 



that gives the corresponding kernel function which defines exactly the same 
search direction and complexity. Nevertheless, for (p (t) = t it corresponds to 
the kernel function k(t) = Ut 1 — 1) — logr of the logarithmic barrier function, 
which gives the classical search direction. This new search direction is also 
analyzed by the authors (Peng, Roos and Terlaky) in a more general context of 
Self-regular proximities measures for linear and (second order) cone optimiza- 
tion problem. See, e.g., [7]. 

The paper is organized as follows. In the next section, the problem is presented. 
Section 3 deals with the new search directions. The description of the algorithm 
and its convergence analysis are presented in section 4. Section 5 ends the paper 
with a conclusion. 

Our notation is the classical one. In particular, R" denotes the space of real 
n -dimensional vectors. Given w, v e R", u T v is their inner product, and || u\\ 
is the Euclidean norm, whereas || uW^ is the /oo-norm. Given a vector u in R", 
U = diag(w) is the n x n diagonal matrix with Uu = w, for all /'. 

2 The statement of the problem 

Through the paper the following assumptions hold. 

• Positive semidefinitness (PSD). The matrix Q is symmetric and positive 
semidefinite. 

• Interior point condition (IPC). There exists (jc°, y°, z°) such that Ax° = 
b, A T y° +z° - Qx° = c, x°> 0, z° > 0. 

• The full rank condition (FR). The matrix A is of rank m (m ^n). 
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The optimal solutions of the primal and the dual are the solutions of nonlinear 
system of (2n + m) of equations: 

Ax = b, x ^ 0 
A T y + z - Qx = c, z ^ 0, (3) 
XiZj = 0, for / = 1, . . . , n. 

The first equation of the system (3) represents the primal feasibility, the second 
one the dual feasibility and the last one the complementarity condition. 

The classical path-following method consists in introducing a positive param- 
eter fi. One considers the nonlinear system parameterized by fx 

Ax = b, 

A T y + z - Qx = c, z > 0, x > 0, (4) 
XiZi = ix, for i* = 1, ... , n. 

It is shown, under our assumptions, that there exists one unique solution 
(x(fx), y(/x), z(a0). The path fx — ► (jc(/x), y(fi), z(fi)) is called the central-path 
(see for instance Wright [8]). It is known that when /x \-+ 0, v(m), z(/x)) 

goes to a solution of (3). Applying the Newton's method for the system (4), 
we develop the classical primal-dual path-following algorithm. 

In the next section, we shall present a new method for approximating the 
central path. 

3 New search directions 

In the proposed method, we replace the n equations jc, y, = n in (4) by n equiv- 
alent equations y ( £ jf L ) = where (p is a real valued function on [0, +00) 
and differentiable on (0, +00) such that (p and (p'(t) > 0, for all / > 0. 

Given > 0 and (x, y, z) such that A T y + z - Qx = c, Ax = b, x > 0 and 
z > 0 but not such that (p (j^j = ^(1) f° r a H '> a new triple (jc + Ajc, y + 
Ay, z + Az) is obtained thanks to the Newton method for solving the nonlinear 
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system. One obtains 

AAx = 0, 

A T Ay + Az- QAx = 0, (5) 

-<p' (—)(ZiAxi +XjAzi) = <p(\)-(p (—] for / = i, — , n. 
fi \ fi J \ M- / 

For commodity let us introduce the vectors v and p of R" defined by 



v, — J and pi = ^-^ — for / — 1 n. 



fx~Zi (p(\)-(p(v-) 

— and pi = j- 

fx Vj(p(vf) 

Next let us introduce the vectors 



d x = X~ l V Ax, d z = Z-'VAz, d y = Ay (6) 
and the matrix 

A = -AV~ ] X and Q = -V~ l XQV~ l X 

then the system reduces to the system 

Ad x = 0, —A T d y — d z + Qd x =0, d x +d z = p, 
which leads to the system 

( v j*f)(i)-{i) 

which is non singular since Q is positive semidefinite and A has a full rank. 
Then d z is computed by the formula 

d z = p -d x . 

Remark 3.1. If we adopt (p(t) = /, then we recover the classical primal-dual 
path-following method 

In this paper, we shall consider <p(t) = <Jt. Then (p'(t) = j^. It follows that 



i=2(l- j^)=2(l-« ( ), 



for all / 

V M ) 
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or 

P = 2{e- v) (7) 
where e = (l,...,l) r and the system (5) becomes 

AAx = 0, 

A T Ay + Az- QAx = 0, (8) 
ZiAxj + XiAzi = fiViPi for /' = 1, . . . , n. 

In addition with the notation in (6), we have 

XiAzi + ZiAxj = fiVi(d Xi +d Zi ) for i = 1, w, (9) 

and 

lid Xi d Zi = Ajc/Az, for i = 1, (10) 

4 Description of the algorithm and its convergence analysis 

Now we define a proximity measure for the central path as 

8(XZe,fi) = ^- = \\e-v\\. (11) 
In addition, we define the vector q by 

q t = d Xi - d Zl . 

Then 

p 2 -q 2 
d xAi = — j — - 

and 

IklKllpll- 02) 

This last inequality follows from 

\\p\\ 2 = \\q\\ 2 +2d T x d z 

since djd z = dJ(Qd x - A T d y ) = d T x Qd x ^ 0, with Ad x = 0 and Q is a 
positive semidefinite matrix. 
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4. 1 The algorithm 

Let e > 0 be the given tolerance, 0 < 9 < 1 the update parameter (default 9 — 
1/(2 y/n)), and 0 < r < 1 the proximity parameter (default x = \). Assume 

0 T 0 

that for the triple (jc°, y°, z°) IPC holds and let fi° = {JL ^ S -- In addition, we 
assume that 8(X°Z°e, ii) < r. 

begin 

x := x°; v := y°; z := z°; 
/x := /x°; 

while x T z = n/ne do 
begin 

fi := (1 

compute (Ax, Ay, Az) from (8) 

x := jc + Ax; 

y := y + Ay; 

z := z + Az; 
end 
end 

In the next section we prove that the algorithm converges quadratically and 
locally and in a polynomial time. 

4.2 Complexity analysis 

In the following lemma, we state a condition which ensures the feasibility of 
the full Newton step. Let .v = v + Aa and z = z + Az be the new iterate after 
a full Newton step. 

Lemma 4.1. Let h — 8(XZe, fx) < 1. Then the full Newton step is strictly 
feasible, hence 

x > 0 and z > 0. 

Proof. For each 0 ^ a ^ 1. Let x(a) = x + a Ax and z(a) = z + a Az. 
Hence 

Xi (a)zi (a) = XjZi + aU, Az, + z, Ax, ) + a 2 Ax, Az, for all / = 1 , . . . , n. 
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Now, in view of (9) and (10) we have 

Xi(a)zi(a) = fi(vf + aviidxi + dzi) + a 2 dxjdzi) for all i = 1, . . . , n. 
In addition from (7), we have 

Vi + ~ = I, for / = 1, . . . , n. 

and thus 

vf + ViPi = 1 — -j- for i = 1, . . 



Thereby 



X|(a)%(«r) = /x((l - a)u f 2 + a (u 2 4- v f p/) + — (p 2 - 4, 2 ) ), 



and 



vf = (1 - a)u f 2 +a/l-(l- a)^ - , for i = 1, . . . , n. 



Thus the inequality x, (a)zi(oc) > 0 holds if 

,2 



max 

i 



Using (11) and (12) we get 

««(|a-«4-4|) < o-«)||ft+HI!lll 

« (1 _ a) I^ +a w 2 



4 



(13) 



^ 4 

Hence, i,-(a)z,-(a) > 0 for each 0 < a < 1. Since Jc, (a) and z, (a) are linear 
function of a, then they do not change sign on the interval [0, 1] and for a = 0 
we have x t (0) > 0 and z,(0) > 0. This leads to i, (l) > 0 and z,(l) > 0 for 
all i. □ 

In the next lemma we state that under the condition 8(XZe, fz) < 1, the full 
Newton step is quadratically convergent. 
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Lemma 4.2. Let x = x + Ax and z = z + Az be the iteration obtained after 
a full step. Suppose 8{XZe, fi) < 1. Then 

8 2 

8{XZe, /J,) ^ 



Thus 8(XZe, fx) < 8 2 , which means quadratic convergence of the full Newton 
step. 

Proof. Setting in ( 1 3) a = 1 , we get 

2 q 2 

vf = l--± for i = (14) 

Hence 
and thereby 



nun 

i 



inv^Jl-8 2 . (15) 



On the other hand, we have 

n 



i=\ 



(i +min / u l -r ~; 



and using (15) we get 



8 2 (XZe, < ! £(1 _ vf) : 
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Using (1 1), (12) and (14) we get 



8'(XZe,fi) ^ 



4 2 (i + vr^y 

(E- =1 ^ 2 ) 2 



Hence 



l / Ml 2 \ 2 

8 2 



Off 



<5(XZe, /x) ^ 



1 + Vl^S 2 "' 

This proves the Lemma. □ 

The next Lemma gives an upper bound for the duality gap after a full Newton 
step. 

Lemma 4.3. Let x = x + Ajc and z = z + Az. Then the duality gap is 

x z = n {n - — J , 



hence 



x T z < \in. 



Proof. We have seen from the previous lemmas that 

v. f = 1 f- for i = 1, . . . , n. 

4 



Hence 



£«-±(i-$)-.-f 



but since 

n 

i r z = ^^5f, (16) 



»=i 

Comp. Appl. Math., Vol. 25, N. 1, 2006 



Copyrighted material 



MOHAMED ACHACHE 



107 



it follows that 



x' Z 



( \\q\\ 2 \ 



This completes the proof. □ 

The next Lemma discusses the influence on the proximity measure of the 
Newton process followed by a step along the central path. 

Lemma 4.4. Let 8 = 8(XZe, /x) < 1 and /x + = (1 -9)fx, where 0 < 9 < 1. 
Then 

9Jh~ + 8 
^ ~ 1 -0 + VO -0)0 -S) 
Furthermore, if 8 < 1/2, 9 = \/(2^/n) and n ^ 4, then we get 

8(XZe, fi+) < i . 

Proof. We have 



8 2 (XZe,fx + ) = 



e - 



1 -9 



x y 



|| Vl -0e-o|| 






1=1 








(Vi 


- e + o f ) 2 

- e + s,) 2 


vM<i-*)-o, 2 ) 2 

£r (vi - <• + so 2 







^((l-9)-v?) 2 
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Hence by using (13) and (14) we get 



8(XZe,n + )K 



Now, in view of (10) and (1 1) we deduce that 

1 



n / 2 \ ' 



8(XZe,fjL + ) ^ 



1 - 6 + y/{\-e)(\ 
1 



1 -(9 + v/(l -6»)(1 -5 2 ) 
Finally, we get 

8(XZe, fi+) ^ . 1 = \0jh~ + S] . 

l-e + y/(\-0)(\-8 2 ) L J 

This completes the proof. □ 

It results from Lemma 4.4 that for the default 0 = l/(2y/n), the (x, z) > 0 
and 8(XZe,fx) < 1/2 are maintained during the algorithm. Hence the algo- 
rithm is well defined. 

In the next Lemma we compute an upper bound for the total number of 
iterations produced by the algorithm. 



Lemma 4.5. Assume that x°and z°are strictly feasible, = and 
8(X°Z°e, n) < j. Moreover, let x k and z k be the vectors obtained after k 

T 

iterations. Then the inequality (x k ) z k ^ s is satisfied for 



k > [e log —r-[ 



Proof. We have after k iterations that fi k = (1 - 0)V°. Using Lemma 4.3, 
it follows that (x k ) T z k < nfi k = (1 - 0)*(jc°) r z o . Hence the inequality 
(x k ) z k ^ £ holds if (1 - 6) k (x°) T z° ^ e. By taking logarithms of both 
sides, we obtain 

*log(l-0) + log(;c o ) r z o <log£. 
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The estimate for the log function 

-log(l -0) >e for all 0 ^ 6 < I. 
Therefore the inequality (x k ) T z* < s is fulfilled if k > | log fi^. □ 

For the default 0 = l/(2 > /»), we obtain the following theorem. 

Theorem 4.6. Suppose that the pair (x°, z°) is strictly feasible and let fM° = 
(x°) T z°/n. If 0 = \/(2y/n), then the algorithm 4.1 requires at most 




iterations. For the resulting vectors we have (xr ) z k ^ e. 
5 Conclusion 

In this paper, we have described a new primal-dual path-following method to 
solve convex quadratic programs. We have proved that the short-update algo- 
rithm can find an e-solution for (QP) in 

2Vnlog — - 

s 

iterations, where x° and z° are the initial positive starting points. The initializa- 
tion of the algorithm can be achieved as in linear optimization by using self-dual 
embedding techniques. Hence, if *° = y° = e, we get the best well-known 
short step complexity 

Iv^log (^j . 
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Abstract. We study the effect of the nonlinear dependence of the iterate x k of Conjugate 
Gradients method (CG) from the data b in the GCV procedure to stop the iterations. We compare 
two versions of using GCV to stop CG. In one version we compute the GCV function with the 
iterate x depending linearly from the data b and the other one depending nonlinearly. We have 
tested the two versions in a large scale problem: positron emission tomography (PET). Our results 
suggest the necessity of considering the nonlinearity for the GCV function to obtain a reasonable 
stopping criterion. 
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1 Introduction 

We consider the problem of estimating a solution x of 



where A is a real m x n matrix, b is an m -vector of observations and € is an 
error vector. We assume that the components of e are random with zero mean, 
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Ax + e = b , 



(1.1) 
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uncorrelated and with variance o 2 (unknown); i.e.; 



Ee = 0, Ee€'=(T 2 I, 



(1.2) 



where E denotes the expectation and / is the identity matrix. 

If A is a full rank matrix, the Gauss-Markov Theorem states that the least 
squares solution of (1.1), i.e., x = (A'A^A'b, is the best unbiased linear 
estimator of x, meaning that it is the minimum variance estimator (see, for 
example [26]). But, if A is ill-conditioned, this minimum variance is still large. 
It is well known that, if we allow the estimator to be biased the variance could 
be drastically reduced (see [28, 29]). 

One way to do this, is by considering solutions of regularized problems of 
the form 



for a positive real number X, where B is an n x n matrix, that introduces a priori 
information on the problem. For example, B could be the matrix associated to 
the discretization of derivatives enforcing smoothness of x. 



Another way is to apply an iterative method that is convergent to the solution 
of 



In this case, a regularization effect is obtained by choosing as 'solution' an 
early iterate of the method (See [31, 9] and the references therein, or see [14, 
Chapter 7]). Here the iteration number k plays the role of the parameter a of 
the first approach. For stationary iterative methods in [10] and [23] it is shown 
that this approach is equivalent in some sense to the first one. 

A crucial point in both approaches is the choice of the regularization parameter 
X in (1.3) or the iteration number k in the second approach. A possibility is to 
approximate a value of a, such that, the average of the mean square error is 
minimum, i.e., X solves the problem 



minimize \\Ax - b\\ 2 + Xx'Bx , 



(1.3) 



minimize ||/4x-fr|| 



(1.4) 



minimize ET(X) , 



(1.5) 



where 



7(a) = -||A;c-A*,|| 2 , 



(1.6) 



m 
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where x k is the solution of (1.3) and x is one solution of (1.1). Probably, the 
most popular of these methods is Generalized Cross-Validation (GCV) (See 
[33, 16]). 

If x k is a estimator for x (the solution of ( 1 .3) or one iterate of a iterative method 
convergent to the solution of (1.4)), the influence operator A x is defined as 

A x b = Ax k . (1.7) 

For the solution of (1.3), the influence operator is given by 

A k b = A(A'A + \B)- ] A'b (1.8) 

If A k is given by (1.8), in [6, 13] the GCV function V{X) was defined as 

i\\b-Ax k \\ 2 
V(A) = -f -j. (1.9) 

GCV chooses the regularization parameter by minimizing V(k). This estimate 
is a rotationally invariant version of Allen's PRESS, or Cross-Validation [1]. 
In [6, 13] was proven that, if A k is given by (1.8), then 

\ET{\)- EV{X) + o 2 \ ( fji 2 \ 



(2m + ^) — ^ (i.io) 

V flij (1 - Ml) 2 



ET(k) 

where /-ij = ^TrA k /x 2 = ^TrAj, whenever 0 < < 1. In [8] Elden presents 
a method to compute the trace-term of the GCV function using bidiagonaliza- 
tion. In [12], Golub and Von Matt proposed an iterative method to approximate 
the trace-term of the GCV function for large scale problems. In [21], Golub, 
Nguyen and Milanfar extended the derivation of the GCV function to under- 
determined systems with applications to superresolution reconstruction. 

Another approach is to apply an iterative method directly to the least square 
problem 

minimize || Ax -b\\ 2 

and take as as a 'solution' an early iterate of the method. In [10] and [23], it is 
proven that this approach is equivalent, in some sense, to ( 1 .3). Now, the role of 
the parameter k is played by the iteration number k and the GCV functional can 
be used to choose this number, that is, as a stopping rule, as shown in [32, 25] for 
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stationary linear iterative methods. It is well known that the Conjugate Gradients 
method applied to the normal equations (CGLS) associated to the problem (1.1) 
can achieve its best accuracy significantly faster than stationary methods and that 
CGLS can be used as a regularizing algorithm (see [22]). However, the iterations 
of CG are a nonlinear function of the data vector b, and, as pointed out in [15], 
a very good stopping rule should be used. In [15], Hanke and Hansen suggest 
a Monte Carlo GCV procedure to stop CGLS and the u-method ([15], pages 
298-299). But the nonlinearity of the influence operator in the case of CGLS is 
not taken into account. While for the v-method the influence operator is linear 
or affine, for CGLS it is nonlinear. However they observe that the application of 
their procedure "can suffer severely from the nonlinearity of CGLS". Also in [ 1 5] 
another even more crude approximation is used. In the later case the denominator 

of the GCV functional is approximated by the expression ^1 — , where 

p is the matrix rank; and this is the same for every value of b ([15], page 307). 

The main goal of this paper is to compare a Monte Carlo GCV procedure 
deduced in Section 2 that takes into account nonlinearity with another one that 
linearizes the influence operator for large scale problems like that arising in 
Positron Emission Tomography (PET). 

Another alternative to GCV can be the use of the L-curve (see [15]), that is 
not considered in this article because of the fact that we wanted to make explicit 
the importance of considering the nonlinear part of CG and because of the lack 
in that approach (the L-curve) of most of the statistical information provided by 
the problem. 

In Section 2 we derive an inequality like (1.10), when A k is a nonlinear 
function of b. We also show how to use a Monte Carlo approach, as in [1 1], 
in order to compute the denominator of our GCV functional. In Section 3 we 
describe the implementation of our procedure applied to a preconditioned Con- 
jugate Gradients Method (PCCGLS). In Section 4 we present numerical experi- 
ments simulating a PET problem and in Section 5 some concluding remarks. 
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2 GCV for a nonlinear influence operator 

We would like to obtain a good estimate for 

T{k) = -\\Ax-Ax k \\\ (2.1) 
m 

using the relative residual defined by 

U(k) = -\\b - Ax x \\ 2 (2.2) 
m 

and our information on the error (1.2). 
The next result is a technical lemma. 



Lemma 1. Let F(k), g{k), G(k), r(k) and H (k) be real valued functions with 
k in R and a real number. If 

G(k) = F(k) + or( 1 - 2g(k)) + r(k) (2.3) 

and 

G(k) 
(1 - g(k))- 
then the following expression is valid: 

"<"-™>-" = ' (- ^^^W lgjk^gjkA (2.5) 

Proof. From (2.3) e (2.4) we have that 

H (a) - F(k) — a 



F(k) 

1 ( F(k)+a(\-2g(k)) + r(k) _ (1 - g(k)) 2 (F(k) + a) \ 
(\-g(k)) 2 \ F(k) F(k) ) 



(2.6) 



that implies (2.5). □ 
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Throughout this paper we will assume that the influence operator, defined by 

A x {b) = Ax kt (2.7) 

is a continuously differentiable operator as a function of b, with b varying in an 
open set £2, that contains Ax, and if we denote by DA k () the Jacobian of A k (-), 
we have the next result. 

Proposition 2. Let {x k } be a family of estimators for the solution of ( 1 . 1 ), and 
(1.2). For each k, A k {-) is continuously differentiable in £2, then 

EU(k) = ET{k) + o 2 (\ - ^Tr(DA k (Ax))j + E&O&U)) , (2.8) 

where U(X) is given by (2.2), T(k) by (2.1) and the function t\(e'e) is such 
that ||0x(e'e)|| < M€ l €, for some M > 0. 

Proof. Expanding the square of the residual of ( 1 . 1 ) and ( 1 .6), at x k we get that 

\\b - Ax k \\ 2 - \\Ax - Ax k \\ 2 = \\Ax + e - Ax k \\ 2 - \\Ax - Ax k \\ 2 

= ||e|| 2 + 2€'(Ax-A*x). 

Now, expanding A x ( ) at Ax + e around the point Ax in Taylor formula: 

A k (Ax + e) = A k (Ax) + DA k (Ax)€ + O^e'e) (2.10) 
where O k {€ x €) satisfies the hypothesis. Using (1.1) and (2.10) above we obtain 
Ax k = A k (b) = A k (Ax + e) = A k (Ax) + DA k (Ax)t + O k (€*€) . (2.1 1) 

Dividing (2.9) by m, applying the expectation and using (1.2), (2.1 1) we obtain 

1 2 

EU(k) = ET(X) + -E(\\€\\ 2 ) + -[E(€ t Ax)-E{6 t Ax k )} 

= ET(k)-^-a 2 --[E(€ t A k (Ax)) 

(2.12D) 

+ E( € 'DA k (Ax)e) + E O k (c'e))] 
= ET(\) + a 2 [\-^Tr(DA k (Ax))^ + E{€'O k (€'6)) , 
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Let we introduce the GCV function for nonlinear influence operators. For 
Ax € Q, we define the GCV functional by 

±\\b- Ax k \\ 2 

V(X) = —-^ (2.13) 

[17>(/ - DA k (Ax))Y 

In spite of the fact that we do not know the value of Ax, we will show bellow 
how we can obtain a good approximation for the denominator of (2.13) for large 
scale problems. 

Theorem 3. Let {x k } be a family of estimators for the solution of ( 1 . 1 ), for 
which (1.2) is valid, and for V(k) given by (2. 13) the following equality holds 

EV(k) - ET(k)-o 2 



ET{\) 

(2.14) 

(1-MiW) 2 V ET(k) + MU / 
w^/x,(X) = £7>(D/U(Ajc))am/r(A.) = E{e< O k (e'e)). 

Proof. By Proposition 2 the inequality (2.8) holds. Applying Lemma 1 we 
get (2.14). □ 

The next Proposition is the basis of our approximation scheme. 

Proposition 4. Let w = (w\, w m ) be a vector of random components 

with normal distribution, that is, w ~ N(0, I inX m)- Suppose that A k (-) is 
continuously differentiable in the open set £2 that contains Ax. Let x\ and x^ the 
estimators corresponding to the data vectors b + 8w and b - Sw respectively. 
Then 



[w-A(x\-x$)/28]\ 



w'w 



E — ; 

(2.15) 

1 „ ^ A ^(w t O k {€ , e + 2&Ww\ + 8 1 w t w) 



= -Tr(I - DA k (Ax)) + E - — ! — ) 

m \ urw / 
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Proof. Expanding A k (-) at Ax + e + Sw and Ax + e + Sw around the point 
Ajc up to the first order: 



Axj 1 = A k (Ax + €+8w) 

= A x (Ajc) + DA k (Ax)(t + + Ox [(6 + «u>)'(€ + Sw)] 

and 

Axj = Ax (Ax + €-8w) 

= A x (Ajc) + DA x (Ajc)(6 - Sw) + O x [(6 - Sw)'(€ - Sw)] . 

Thus, we obtain 



(2.16) 



(2.17) 



4 Ax^~ 

— !— — 2 - = DA x (Ax)w + O k (c'c + 25|€ , u>| + <5 2 w;'u;) (2.18) 
2d 

and the result follows using a slight extension of the Theorem 2.2 of [1 1]. □ 

For the sake of consistency with the usual mathematical notation for the it- 
erations of an algorithm, in what follows, the parameter X will be substituted 
by k. 

3 Monte Carlo Implementation of GCV for Conjugate Gradients 

Our main goal in this section is to apply our Monte Carlo Implementation of 
GCV as stopping rule for Conjugate Gradients method for solving (1.1). Now, 
each iterate x k is an estimate of the solution of (1.1) and k plays the role of A. 
In Section 2 we showed that GCV is applicable to influence operators that are 
continuously differentiable in an open set that contains Ax. 

The Conjugate Gradients method of Hestenes and Stiefel [3] applied to the 
normal equations (CGLS) can be written as: (for a given x°) 

r° = b-Ax°, s° = A f r°, w° = s° and for fc = 0,l,... 
p k = Aw k 
k \\s k \\ 2 



a = 



*ll2 



x k+l _ x k _(_ Qf* u^* 
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r k+i _ r k _ a k pk 
g k+i _ A t r k+] 



\\s k \\ 2 
w k+] =s k+] +0 k w k 

For each k > 0, the operator T(k)(b) = x k is nonlinear and it is continuously 
differentiable outside the termination set (closed) 

M k = {fc € R m |P:R(y4)b - Ax° is a linear combination of k eigenvectors of A A 1 J , 

where Py denotes the orthogonal projection onto the subspace U (i.e., £2 = 
(R n - M k )). If A is an ill-conditioned operator, then T(k)(-) is discontinuous in 
, as pointed out by Eicke, Louis and Plato in [7] in the infinite dimensional 
context. The same arguments of [7] can be easily used to understand the behavior 
of the method in the discrete case. For b in A/*_|, we define b l = b + o\*Uu 
where u\ is the singular vector associated to the singular value a/. From the fact 
that b l e M k it follows that 

\\T(k)(b) - T(k)(b l )\\ = \\A f b - A f b'\\ = -V. 

°C 

Thus, if or/ ^ 0, then \\T(k)(b) - T(k)(b l )\\ * oo. And if b e M k ^ and/?,,, is/? 
perturbed by a random vector w, we should also have \ \T(k)(b) — T( k )(b w ) \\ ~ 
oo. But, usually in ill-posed problems the vector b has all the components of 
the singular vectors. It means that b does not belong to a termination set M k , 
for any k. 

Taking into account the results of the preceding section we describe the algo- 
rithm for calculating an approximation of the GCV functional (2.13) by using 
central differences to approximate the direcional derivative. 

Algorithm 1. 

(i) Generate a pseudo-random vector w = (w\, . . . , w m )' e R m from a nor- 
mal distribution with standard deviation equal to one, i.e., 
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(ii) For each k, compute the iterates x k x k and x\ corresponding to b, b — Sw 
eb + 8w, respectively, where 8 = 10~ 4 . Take 

™>-( " l "- A( i-* m )\ 0.19) 

as an approximation for (^Tr [I m - DAk{Ax)])" ; 

(iii) For each k take 

M\b- Ax k \\ 2 
V ^ k) = \ /m < 3 - 20 > 

<t>NL(k) 

4 Positron emission tomography 

The goal of positron emission tomography (PET) is the quantitative determina- 
tion of the moment-to-moment changes in the chemistry and flow physiology 
of radioactive labelled components inside the body. The mathematical problem 
consists of reconstructing a function representing the distribution of radioactivity 
in a body cross-section from measured data that are the total activity along lines 
of known location. One of the main differences between this problem and that 
arising in X-ray tomography [17] is that here measurements tend to be much 
more noisy, so, direct inversion using convolution backprojection (CBP) doesn't 
necessarily give the best results (see [30]). 

In positron emission tomography (PET) [27], the isotope used emits positrons 
which annihilate with nearly electrons generating two photons travelling away 
from each other in (nearly) opposite directions; the number of such photons pairs 
(detected in time coincidence) for each line or pair of detectors is related to the 
integral of the concentration of isotope along the line. 

Suppose now that we discretize the problem by subdividing the reconstruction 
region into n small square-shaped picture elements (pixels, for short) and we 
assume that the activity in each pixel / is a constant, denoted by xj. If we 
count bj coincidences along m lines and a,, denotes the probability that a photon 
emitted by pixel j is detected by pair i, then v, is a sample from a Poisson 
distribution whose expected value is ]T^=i OijXj. 

It is well known that ART, with small relaxation parameters [18], approximates 
a weighted least squares solution of ( 1 . 1 ). This fact suggests that the application 
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of Conjugate Gradients to the system (4.21) could give similar or better results 
[24] and it is a reasonable example to test the capability of GCV as a stopping 
rule for CG. 

With the objetive of comparing the GCV procedures in Computerized To- 
mography we apply the Conjugate Gradients method "preconditioned" with a 
symmetric ART (from algebraic reconstruction technique) method, as presented 
recently in [24], i.e. we apply CG to solve the generalized least squares problem 
(CGGLS) 

minHC-'Ob- Ax)\\ 2 . 

where = (D + (oL)D~ s, if A A 1 = L + D + U, where L is the strictly lower 
triangular part and D the diagonal of A A' respectively. 
This corresponds to apply CG to the system 

A'CZ'CZ l Ax = A t CZ t CZ l b, (4.1) 

In the ART method, co is the relaxation parameter. Here it has a different role, 
it is the weight of the lower triangular part of A in the preconditioning. When co 
is zero, this corresponds to a diagonal scaling. 

The Conjugate Gradients method applied to (4.21) can be written as: 

r° = C~\b - Ax°), s° = A'C-'r\ w° = s° and for k = 0, 1, . .. 
P k = C~'Aw k 

IIP*II 2 

x k+\ _ x k + a k w k 
*+l k k k 

r s= r — a p 
s k+l = A'C-'r k+l 

, = ll^ +1 ll 2 

ll^ll 2 
w k+i _ s k+\ + pk w k 

The elements of the matrix L are not explicitly known, but in [2] was shown 
an efficient way to compute the products A'C'Jr and C~ x Aw. 
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In our numerical experiments we used the programming system SNARK93, 
developed by the Medical Image Processing Group of the University of Penn- 
sylvania [5]. The images to be reconstructed (phantom) were obtained from a 
computerized atlas based on typical values inside the brain, as in [18]. The data 
collection geometry was a divergent one simulating a typical PET data acqui- 
sition [5]. We used a discretization with n = 95 x 95 pixels and the divergent 
geometry had 300 views, of 101 rays each, a total number of m = 30292 equa- 
tions (8 rays were not considered because their lack of intersection with the 
image region). The starting point was a uniform image x° — (a, . . . , a)', where 
a is an approximation of the average density of the phantom given by 

a = x^m v^n ' (4-2) 

The choice of a uniform non-zero starting point was advocated by L. Kaufman 
[20] and it is widely accepted as the best choice for many researchers in PET. 
The vector b was taken from a pseudorandom number generator with a Poisson 
distribution (see [ 1 9]). The total photon count was 2 022 085, 99 1 1 79, 5 14 925 
and 238 172 (monotonically increasing noise). 

We have performed experiments applying Algorithm 1 with photon counts 
2 022 085, 991 179, 514 925 and 238 172 and with six values of co, between 0.0 
and 0.025. For each value co and each photon count we repeated the experiments 
ten times. The ten tests gave very similar results. 

We plotted, for one of the tests, in Figure 1 the functions 

lOO^H** -*|| 2 , ^±^\\Ax k - Ax\\\ Vmc-100 and V HH - 100 

against the iteration number k, for co = 0.0 (diagonal scaling) and a> = 0.025. 
Here we call Vhh the Monte Carlo GCV function, if we suppose that PCCGLS 
is a linear function of the data vector b, as proposed by Hanke-Hansen in [15]. 
As expected, the minimum of our Monte Carlo GCV function Vmc(^) coincides 
with that of 1 1 A (x k — x) | | 2 , and also the curves are very similar. We can see that 
using Vhh can give wrong results in the case of Conjugate Gradients. Figures 
2 and 3 show the reconstructions corresponding to the curves in row 1 of the 
Fig. 1. 
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iteration, k iteration, k 

Figure 1 - Loss functions for the method CGGLS with co = 0.0 (diagonal scaling) 
(left) and co = 0.025 (right). The rows 1-4 correspond to 2 022 085, 991 179, 
514 925 and 238 172 photon counts, respectively. The solid thin line. 1 correspond 
to IOO^jIIjc* — x\\ 2 ; the thick ones, to 2^92 \\Ax k — j4jc|| 2 ; the dashed lines to transla- 
tions of Vmc and dotted ones to translations of Vhh, if we suppose that CGGLS is 
a linear function of b. 
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Figure 2 - A 95 x 95 phantom randomly generated and its reconstructions using CGGLS 
with a> = 0.0 (diagonal scaling). 
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Figure 3 - A 95 x 95 phantom randomly generated and its reconstructions using CGGLS 
with (o = 0.025. 
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5 Concluding remarks 

Iterative methods in Emission Computed Tomography, like RAMLA [4], that 
are being used nowadays by PET scanners (see for example http://www. 
medical.philips.com/us/products/pet/products/cpet/), are fast and produce high 
quality pictures (because they take into account the Poisson nature of the noise) 
in few iterations. So it was not our goal in this paper to compare Conjugate 
Gradient with those methods, but to show that our approximation of GCV is 
applicable as a stopping rule to the Conjugate Gradient method for large scale 
ill-posed problems, similar to Positron Emission Tomography (very large scale, 
not severely ill-posed and relatively large noise in the data). Some authors (see, 
for example, [15, 9]) have suggested the use of other approximations to GCV. 
Our experimental results show (Fig. 1) that Algorithm 1 gives very good (much 
better than Hanke et al's) approximations for the minimum of \\A(x k - x)\\ 2 . 

Further research is needed to extend these applications of GCV to more gen- 
eral nonlinear methods and nonlinear problems. We also need more specific 
asymptotic theoretical results. 
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